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ABSTRACT 


This  thesis  considers  the  problem  of  time-optimal  spacecraft  slew 
maneuvers.  Since  the  work  of  Bilimoria  and  Wie  it  has  been  known  that  the  time- 
optimal  reorientation  of  a  symmetric  rigid  body  was  not  the  eigenaxis  maneuver 
once  thought  to  be  correct.  Here,  this  concept  is  extended  to  axisymmetric  and 
asymmetric  rigid  body  reorientations  with  idealized  independent  torque 
generating  devices.  The  premise  that  the  time-optimal  maneuver  is  not,  in 
general,  an  eigenaxis  maneuver,  is  shown  to  hold  for  all  spacecraft 
configurations.  The  methodology  is  then  extended  to  include  spacecraft  control 
systems  employing  magnetic  torque  rods,  a  combination  of  pitch  bias  wheel  with 
magnetic  torque  rods,  and  finally  to  control  systems  employing  single  gimbal 
control  moment  gyros.  The  resulting  control  solutions,  designed  within  the 
limitations  of  the  actuators,  eliminate  the  requirement  to  avoid  actuator 
singularities.  Finally,  by  employing  sampled-state  feedback  the  viability  of  real¬ 
time  optimal  closed  loop  control  is  demonstrated. 
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I.  INTRODUCTION 


A.  MOTIVATION 

Modern  satellites  require  more  precision  and  agility  than  ever  before.  In 
both  military  and  commercial  applications  the  ability  to  rapidly  maneuver 
represents  an  increase  in  mission  effectiveness  and  productivity.  Rapid 
retargeting  maneuvers  translate  directly  into  more  time  on  the  intended  object 
and  more  observations  per  orbit.  With  commercial  Earth-imaging  satellites  like 
Ikonos  and  research  into  agile  microsatellites1  on  the  rise,  the  need  for  time- 
optimal  satellite  control  is  greater  than  ever.  Planned  future  satellite  missions  in 
support  of  missile-defense  and  Earth  observation  will  rely  on  satellite  agility  and 
minimum-time  maneuvers  for  mission  success. 

B.  THE  PROBLEM 

Spacecraft  time-optimal  attitude  maneuvers  have  held  the  interest  of 
engineers  and  mathematicians  for  decades.  In  their  paper,  “Survey  of  Time- 
Optimal  Attitude  Maneuvers,”  Scriverner  and  Thompson  provide  a  summary  of 
work  in  this  and  related  areas2.  They  also  present  the  principal  difficulties 
associated  with  solving  the  time-optimal  attitude  maneuver  problem. 

Unlike  the  pointing  problem  which  lends  itself  well  to  linearization 
techniques,  slew  problems,  especially  large  angle  slew  maneuvers  are  highly 
nonlinear.  Euler’s  equations,  which  represent  the  rotational  motion  of  a  rigid 
body,  consist  of  three,  coupled,  nonlinear  differential  equations. 

l/0X+(l2-ly)(0/0y=Ui 

/y®y  +  (  4  -  4  )®/°z  =  U2 

/zC6z  +  (/y-/x)cOyCOx=t73 

Additionally,  the  high  angular  velocities  associated  with  time-optimal 
maneuvering  cause  significant  gyroscopic  stiffness  through  the  non-linear  terms. 


1 


While  these  equations  completely  describe  the  rotational  motion  of  the 
body  with  respect  to  an  inertial  frame  they  do  not  determine  the  spacecraft’s 
attitude.  To  describe  the  attitude  of  the  spacecraft  Euler  parameters  are 
generally  used  since  they  provide  singularity  free  kinematics: 
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These  nomlinear  ordinary  differential  equations  have  no  closed-form  solutions 
except  for  a  small  number  of  cases  involving  simple  rotations  or  simple  geometry 
(i.e.,  torque  free  motion  of  an  axisymmetric  body3  4). 

The  time-optimal  nature  of  the  problem  adds  additional  difficulties.  Major 
problems  are  encountered  in  solving  the  boundary  value  problem  that  arises  from 
the  application  of  the  Minimum  Principle.  The  problem,  formulated  in  Chapter  II, 
has  no  known  numerical  or  analytical  solution  except  when  certain  restrictions 
are  applied.  These  assumptions  have  taken  the  form  of  specific  configurations  or 
restricted  motions5. 


C.  HISTORICAL  BACKGROUND 

A  short  summary  of  the  historical  work  that  influenced  this  research  is 
presented  here.  Beginning  with  their  landmark  work  in  1993,  Bilimoria  and  Wie 
demonstrated  with  extensive  analytical  modeling  and  numerical  analysis  that  the 
time-optimal  maneuver  was  not  the  long-assumed  eigenaxis  maneuver3.  Their 
methods  and  results  are  examined  in  detail  in  Chapter  III.  Bilimoria  and  Wie  later 
extended  this  work  to  an  axisymmetric  body  and  observed  the  same 
characteristic  precession  they  had  noted  earlier.  Beyers  and  Vadali8  reproduced 
the  results  of  Bilimoria  and  Wie  but  focused  on  developing  a  control  algorithm. 
They  used  linearization  and  the  switch  time  optimization  (STO)  algorithm 
developed  by  Meier  and  Bryson9  to  produce  an  algorithm  that  could  be 
implemented  in  real-time. 
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In  the  context  of  magnetic  control,  Junkins  et  al.10  modeled  a  single 
controller  aligned  with  the  spin  axis  of  a  spin-stabilized  symmetric  spacecraft. 
The  time-optimal  solutions  were  found  through  an  interactive  graphical 
technique.  These  results  were  eventually  implemented  for  open  loop  control  on 
the  NOVA-1  spacecraft11. 


D.  CONCLUSION 

Solving  the  problem  of  spacecraft  time-optimal  control  has  occupied  the 
interest  of  engineers  and  mathematicians  for  years.  The  problem  is  simple  to 
formulate  and  yet  solutions  have  been  difficult  to  obtain.  In  this  thesis  we  extend 
the  present  body  of  work  to  all  spacecraft  moment  of  inertia  configurations.  The 
ideal  actuators  often  studied  to  simplify  the  problem  formulation  are  replaced  with 
magnetic  torque  rods  and  a  combination  of  toque  rods  and  a  pitch-bias  wheel. 
Finally,  the  control  moment  gyro  configuration  most  studied  for  its  singularity 
problems  is  examined  and  shown  to  be  singularity  free  in  the  time-optimal  result. 
This  unpredicted  benefit  of  the  formulation  has  the  potential  to  make  future  work 
in  actuator  singularity  avoidance  moot.  Finally,  computational  speeds  are  shown 
to  be  such  so  as  to  allow  time-optimal  solutions  to  the  nonlinear  plant  to  be 
generated  in  real-time. 
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II.  OPTIMAL  CONTROL  PROBLEM  FORMULATION 


A.  INTRODUCTION 

In  this  section  the  optimal  control  problem  formulation  is  developed.  This 
formulation  will  be  referred  to  throughout  the  remainder  of  this  work.  There  are 
numerous  excellent  books  and  papers  on  optimal  control  theory  which  explain 
these  concepts  in  much  greater  detail.  Interested  readers  are  referred  to  the 
references  [1 , 2,  3,  &  4]  for  further  information. 


B.  METHODOLOGY 

Over  the  years  many  different  methods  of  solving  optimal  control  problems 
have  been  developed.  These  are  broadly  grouped  into  two  categories:  indirect 
and  direct  methods.  In  indirect  methods,  the  necessary  conditions  for  optimality 
are  derived  from  Pontryagin’s  Minimum  Principle  and  solved  to  obtain  the  optimal 
trajectory.  These  methods  are  notoriously  labor  intensive.  In  direct  methods,  the 
optimal  control  problem  is  discretized  into  a  parameter  optimization  problem. 
The  resulting  nonlinear  programming  problem  can  then  be  solved  by  standard 
nonlinear  programming  means. 

In  this  work  we  will  employ  a  Legendre  Pseudospectral  method  encoded 
in  the  reusable  software  package  DIDO5.  Pseudospectral  methods  are  well 
known  in  the  field  of  fluid  dynamics  where  they  are  used  to  numerically  solve 
partial  differential  equations.  Unlike  other  methods  which  employ  piecewise- 
continuous  polynomials  Pseudospectral  methods  are  unique  in  their  application 
of  global  orthogonal  polynomials  as  trial  functions. 

In  the  Legendre  Pseudospectral  method  the  time  domain  of  the  problem  is 
discretized  at  a  special  set  of  Legendre-Gauss-Lobatto  (LGL)  points.  Polynomial 
approximations  of  the  state  and  control  variables  are  considered  where  Lagrange 
polynomials  are  the  trial  functions  and  the  unknown  coefficients  are  the  values  of 
the  state  and  control  variables  at  the  LGL  points  (nodes).  Using  the  properties  of 
the  Lagrange  polynomials  the  nonlinear  differential  state  equations  are 
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transformed  to  nonlinear  algebraic  equations.  These  equations  are  then  posed 
as  a  nonlinear  programming  problem  and  a  sparse  numerical  optimizer  is  used  to 
solve  the  problem.  In  addition,  a  relationship  between  the  costate  variables  and 
the  Lagrange  multipliers  called  the  Covector  Mapping  Theorem,  allows  for 
determination  of  the  costates  at  the  LGL  points6.  This  provides  numerical 
information  that  is  used  to  validate  the  solution’s  optimality.  Interested  readers 
should  refer  to  references  7,  8,  9,  and  10  for  details  regarding  Pseudospectral 
Methods. 

Numerical  results  throughout  this  work  are  specified  in  terms  of  the 
number  of  LGL  points  used  to  obtain  the  solution.  Initial  solutions  were  generally 
based  on  30  LGL  points  with  initial  guesses  restricted  to  two  vectors  representing 
the  initial  and  final  conditions.  Initial  control  solution  guesses  were  arbitrary. 
Where  greater  accuracy  was  desired  the  30  LGL  point  state  and  control  solution 
was  used  as  a  guess  for  a  second  solution  based  on  1 00  LGL  points. 

C.  OPTIMAL  CONTROL  PROB LEM  FORMULATION 

In  an  attempt  to  develop  the  notation  and  methodology  that  will  be  used 
throughout  this  work  we  consider  the  following  optimal  control  problem. 
Determine  the  control  function  u*  (t)  and  the  corresponding  state  trajectory 

x*  (i f )  that  minimize  the  Bolza  cost  functional* 


JU(),u()Jl)  =  E{x(tl),tl)  +  ]F(xl,t),u(t),t)dt 

1 0 

where  xe  Rn  and  ueRm  are  subject  to  the  differential  constraint 

x=f{x{t),Ut))  te[t0,tf] 

and  the  boundary  conditions 

§o{x{t0),t0)=0  e0eRp  and  p<n  +  1 
§f{x{tf),tf)=  0  efeRq  and  q<n  + 1 

(•) :  Functional  dependence  not  specified. 
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and  control  inequality  constraints  of  the  form 

h{u,t)<0  heRr 

where  dh/du  has  full  rank. 

The  Lagrange  multiplier  theory  allows  us  to  adjoin  the  state  equations  and 
constraints  to  the  cost  functional  to  form  an  augmented  cost  functional  as  follows 

J  =  E  ( x{  tf ) ,  tf )  +vj e0  (( x(t0 ) ,  t0 )+  vfT e,  (( x{ tf ) ,  tf )  ■ + 

|  { F(x,  u,t)  +  X{t)T{f<^(,u,t)-x)  +  [i{t)Th{u,t)}dt 

h 

Defining  the  Hamiltonian  as, 

H(X,x,u,t)  =  F(x,u,t)  +  XTf(x,u,t) 

Pontryagin’s  Minimum  Principle  provides  the  following  necessary  conditions  for 
u*  to  be  an  optimal  control. 

H(V,x*,u*,t)  <H(V,x*,u,  t)  ->Hamiltionian  Minimization 

X*  =  ->  Adjoint  equations 

dx 

do 

X(t0)  =  -vT — ?-->  Initial  transversality 

dx0 

X(tf)  =  ^-+vT^L->  Terminal  transversality 
dxf  oXf 

H\tf]  +  —+vT  —  =  0  ->  Hamiltonian  Value 
L  dt  dt 

For  the  case  in  which  the  cost  function  and  constraints  are  linear  in  the 
state  and  control  variables,  no  minimum  exists  for  the  Hamiltonian  minimization 
unless  inequality  constraints  are  imposed  on  the  state  and/or  control  variables.11 
In  the  case,  where  the  inequality  constraints  are  linear  and  placed  only  on  the 
control  variables  the  problem  is  a  special  case  of  linear  programming  problem 
covered  in  detail  by  Bryson^  and  others13.  The  important  result  is  the  application 
of  the  Karush-Kuhn-Tucker  (KKT)  Theorem  and  complementarity  conditions. 

Then,  if  a  minimum  exists,  it  will  require  the  control  variable  to  be  at  a  boundary 
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of  the  feasible  control  region.  Since  the  Hamiltonian  is  subject  to  an  inequality 
constraint  on  the  control  variable  we  apply  the  Karush-Kuhn-Tucker  Theorem 
and  form  the  Lagrangian  of  the  Hamiltonian: 

H  =  H+\xTh  (2.1) 


where  p,  is  a  KKT  multiplier.  Then  the  KKT  Theorem  gives  the  following 
necessary  conditions  for  the  optimal  trajectory 


dH__dH_ 

du,  du 


+ 


a/? 

du 


(1  =  0 


(2.2) 


In  addition,  the  multiplier-constraint  pair  must  satisfy  the  complementarity 
conditions  of  the  KKT  Theorem  which  states: 


<0  hi(u,t)  =  h^(t) 

=  0  if  /7f(0</7,(M</f(0 
>0  hi(u,t)  =  h^J(t) 

unrestricted  h-(t)  =  h“(t) 


(2.3) 


In  subsequent  chapters,  once  the  state  dynamics  of  the  problem  have 
been  established,  be  optimal  control  problem  will  be  presented  in  this  format. 
The  candidate  optimal  control  solution  will  then  be  propagated  through  the  state 
dynamics  to  verify  the  feasibility  of  the  solution.  Verifying  that  this  candidate 
solution  is  optimal  is  more  difficult.  Recall  that  the  Minimum  Principle  supplies 
necessary  conditions  not  sufficient  conditions  for  optimality.  However,  since  the 
Legendre  Pseudospectral  Method,  through  the  covector  mapping  theorem, 
provides  costate  information  at  the  nodes  we  can  evaluate  the  necessary 
conditions  for  optimality.  Where  available,  results  will  be  compared  to  published 
works.  Finally,  engineering  and  physical  insight  will  be  used  to  establish  that  the 
solutions  obtained  are  optimal. 
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III.  IDEALIZIED  TORQUE  ACTUATOR  CONTROL  PROBLEM 


A.  INTRODUCTION 

In  this  section,  we  develop  the  time-optimal  rotational  maneuver  for  a 
rigid-body  with  ideal,  bounded  actuators.  The  idealized  actuator  will  be  defined 
by: 

Torqueout  =[l]u 

where  \u\  <  1 ,  is  the  control  vector  and  /  is  the  identity  matrix.  In  the  case  of  the 

ideal  actuator,  the  control  vector  is  equal  to  the  output  torque  vector.  This 
distinction  allows  for  later  definition  of  the  control  vector  u  as,  for  example, 
magnetic  dipole  moment  (See  Chapter  IV).  In  this  case  the  actuator  is  not 
considered  ideal.  The  closest  physical  approximation  to  the  idealized  actuator  is 
a  thruster. 

Bilimoria  and  Wie  showed  that  the  time-optimal  solution  for  the 
reorientation  of  an  inertial  symmetric  body  was  not  necessarily  an  eigenaxis 
maneuver.i  Shen  and  Tsiotras  examined  time-optimal  reorientations  of 
axisymmetric  spacecraft  using  two  controls.2  Finally,  Proulx  and  Ross  examined 
the  control  structure  and  evaluated  time-optimal  reorientations  of  asymmetric 
spacecraft3.  We  begin  with  a  reexamination  of  the  inertial  symmetric 
reorientation  problem  in  order  to  establish  the  methodology.  Then  the  principles 
are  extended  to  axisymmetric  and  finally  to  asymmetric  spacecraft  in  the  orbital 
frame. 

B.  INERTIALLY  SYMMETRIC  RIGID-BODIES 

Inertial  symmetric  rigid-body  reorientations  represent  the  simplest  problem 
that  we  can  pose.  Yet,  prior  to  the  work  of  Bilimoria  and  Wie  the  solution  was 
misunderstood.  The  eigenaxis  maneuver  about  the  control  axis  was  thought  to 
be  the  time-optimal  maneuver  for  symmetric  bodies  and  near  time-optimal  for 
other  configurations.4  Through  analytical  analysis  and  numerical  simulation 
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Bilimoria  and  Wie  demonstrated  that  this  was  not  true  for  symmetric  bodies. 
They  showed  that  ail  three  control  components  cannot  be  simultaneously  zero  on 
the  time-optimal  trajectory. 


1.  Notation  and  Transformations 

In  order  to  progress  with  the  problem  formulation  we  establish  the 
standard  notation  and  rotation  sequences  used  in  the  development.  The 
spacecraft  will  have  an  assumed  standard  orbit  frame  defined  as: 

Z0  —  Nadir  pointing 

X0  -  Velocity  vector 

Y0  -  completes  the  right  hand  set 

We  choose  the  rotation  sequence  for  the  body  to  orbit  frame 
transformation,  represented  as  Euler  angles  as  \p  -^>0  -Mp  with  axes  order  of 
rotation  3  -^2  ->1 .  Therefore,  \\t  is  the  first  angular  rotation  about  the  zbody 

axis.  The  second  rotation  is  about  the  once  displaced  ybody  axis  by  an  angle 
0  .  The  final  angular  rotation  cp  is  about  the  x-body  axis. 

Then,  following  the  notation  convention  of  Kane5  the  following  angular 
velocities  are  defined: 


N-b 

COg 


— >  Angular  rate  of  body  with  respect  to  inertial  in  body  frame.  (3.1 ) 


N—O  A 


CO 


o 


0 

-®o 

0 


-» Angular  rate  of  orbit  with  respect  to  inertial  in  orbit  frame.  (3.2) 


o—b 

CDo 


co1 

C02 

®3 


-» Angular  rate  of  body  with  respect  to  orbit  in  body  frame.  (3.3) 
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and: 


N65bB  =  %YB°  +  %JB  (3.4) 

The  rotational  transformation  from  the  orbit  frame  to  the  body  frame  is 
referred  to  as  the  direction  cosine  matrix  (DCM).  It  is  defined  from  the  rotation 
sequence  above  and  represented  in  terms  of  the  quaternion  vector  as,6 

4  -  4  -  4  +  4  2 :{q,qz  -  q3q4)  2{q,q3  -  q2q4) 

BC°=  2 (q,q2-q3q4)  q\  -  q\  -  q23  +  q\  2{q2q3-qq4)  (3.5) 

2(q1<73  - q2q4)  2 (q2q3-q,q4)  q\-q\-q\  +q24_ 

Therefore,  the  angular  velocity  of  the  body  with  respect  to  the  Newtonian  frame 
may  be  written  as: 

XJ  =  °«B  +  BC°  %°0  (3.6) 

Then  substituting  from  equation  (3.2)  we  have, 

N6JbB  =  °65bB  -co0Ci2  (3.7) 

where  co0  is  the  magnitude  of  the  orbital  angular  velocity*  and  C/2 is  the  second 

column  of  the  DCM.  Then  by  substituting  the  values  established  in  equation 
(3.5)  we  obtain  the  following  relationships, 

“xl  fal  r  2 iq,q2+q3q4)  “ 

wy  =  co2  -co0  ql~q]  -qz3+  q\  (3.8) 

_©zj  L®3  J  L  _ 


Orbital  angular  velocity  coo  is  sometimes  referred  to  as  “mean  motion”  for  circular  or 
elliptical  orbits. 
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and, 


co1 

V 

"  2(g1<72  +q3q4)  “ 

co2 

= 

®y 

+  co0 

2  2  2  2 

q2  -  -  Qs  +  Qa 

®3_ 

_  2(q2q3  +  QiQ4)  . 

2.  Problem  Formulation 

The  inertial  symmetric  body  is  shown  in  Figure  1 .  The  state  of  the  system 
can  be  completely  defined  by  its  attitude  and  the  time-rate -of-change  of  attitude. 


Figure  1  Inertial  Symmetric  Body 

For  mathematical  simplicity  in  extension  to  future  work,  the  attitude  is 
represented  as  a  quaternion.  Quaternions  have  none  of  the  inherent  singularities 
that  are  well  known  in  other  representations  of  spacecraft  attitude.  Additionally, 
they  are  computationally  more  suited  for  on-board  real-time  processing  since 
there  are  no  trigonometric  functions  to  be  evaluated  in  the  quaternion  kinematics 
equations. 

Euler’s  rotation  theorem  states  that  a  rigid  body  can  be  changed  from  any 
given  initial  orientation  to  an  arbitrary  final  orientation  by  a  single  rotation  about 
an  axis  that  is  fixed  to  the  body  and  stationary  in  an  inertial  reference  frame.7 
This  axis  which  remains  unchanged  in  both  the  body  and  reference  frames  is 
called  the  eigenaxis  or  Euler  Axis.8 
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Defining  the  eigenaxis  in  the  body  frame  as: 

e  =  ej  +  ej  +  e3k 

allows  us  to  define  the  quaternion  vector  as  follows: 

^  =  e.sin(^) 
q2=e2  sin(^) 

^3  =  e3sin  {%) 

Qa  =cos[%) 

where  <|)  is  the  rotation  angle  about  the  eigenaxis.  The  state  of  the  spacecraft  is 
then  represented  by:t 


to 


The  state  dynamics  include  both  quaternion  kinematics  and  rotational 
dynamics.  The  quaternion  kinematics  are  well  known  and  are  repeated  here  for 
completeness. 


(3.10) 


(3.11) 


q  =  with,  £1  = 


0 

CO 

3 

-co2 

co1 

-co3 

0 

0^ 

CO, 

c 

co2 

— 

0 

CO, 

— co1 

— co2 

-co3 

0 

(3.12) 


Termwise,  we  have, 


t  We  have  adopted  the  convention  that  q4  is  the  scalar  quantity  of  the  quaternion  vector. 
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4  =^[®iQ4-«bQ3+c%Q2] 
4  =^KQ3+C^C?4-“3Qi] 

4  =^[-®iQ2  +c°2Qi  +CO3Q4] 
4=|[-C0i4-“24-“34] 


(3.13) 


In  order  to  maintain  consistent  notation  throughout  this  work  we  have  used  the 
subscripts  indicating  that  the  angular  velocities  are  of  the  body  with  respect  to  the 
orbit  in  the  body  frame.  For  this  example,  we  consider  the  spacecraft  located  in 
inertial  space.  Mathematically,  orbital  angular  velocity  is  zero  and  we  have, 


CO. 

CO 

1 

X 

co2 

= 

°>y 

®3_ 

This  distinction  will  become  important  for  the  later  case  where  the  spacecraft  is  in 
orbit  about  the  Earth  centered  inertial  frame. 

The  rotational  dynamics  of  a  rigid  body  are  obtained  by  equating  the 
applied  torque  about  the  center  of  mass  to  the  time  rate  of  change  of  the  angular 
momentum.9  This  well-known  fact  is  expressed  as: 


where  H  is  the  angular  momentum  vector  of  the  rigid-body  about  its  center  of 
mass  with  respect  to  an  inertial  frame  and  M  is  the  external  moment  acting  on 
the  body  about  its  center  of  mass. 

Since  this  is  the  logical  extension  of  Newton’s  2nd  law  it  follows  that  the 
time-derivative  of  angular  momentum  must  be  with  respect  to  an  inertial  frame. 
The  angular  momentum  is  the  product  of  the  moment  of  inertia  and  the  angular 
velocity, 

H  =  ia 3 
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This  allows  us  to  express  Euler’s  rotational  equations  of  motion  for  a  rigid-body  in 
matrix  notation. 

/co  +  cox  /co  =  Mext 

If  we  express  the  moment  of  inertia  and  angular  velocity  in  the  principal  axis 
frame  as: 

"/*  o  o“ 

1=  0  ly  0 

o  o  /z 

co  =co1j61  +co>4  +  g>34 

Then  Euler’s  equations  can  be  expanded  to: 

=/,0)x+(/z-/y)Wy0)z 

M2  =/ycoy  +  (/x-/z)coxcoz  (3.14) 

/W3  =//6z  +  (/y  -/x) wxcoy 

These  equations  are  well  known  and  generally  referred  to  as  Euler’s  Moment 
Equations.10  For  the  case  of  a  symmetric  body  the  gyroscopic  terms  in  equations 
(3.14)  are  zero  and  the  rotational  dynamics  can  be  written  as: 

.  M  .  M  .  M. 

C°y=-T  &z=-r 

x  y  z 

This  reduces  Euler’s  Moment  Equations  from  three-coupled  non-linear  ordinary 
differential  equations  to  three  uncoupled  linear  ordinary  differential  equations. 
Then  normalizing  the  moment  by  the  inertia,  without  any  loss  of  generality  we 
have, 

c bx  =  c6v  =  M2  o T  =  M3  (3.15) 

Taken  together,  equations  (3.13)  and  (3.15)  form  the  system  dynamic 
constraints.  The  normalized  external  torque, Mn  is  chosen  as  the  control 
parameter  and  it  is  limited  somewhat  arbitrarily  to  unity. 
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2.  Time  Optimal  Maneuvers 

With  the  dynamic  constraints  defined  and  presuming  initial  conditions  and 
desired  final  conditions  are  known,  we  can  formally  state  the  time-optimal  control 
problem  for  the  rest-to-rest  maneuver  as  follows: 

Minimize  J(x{),u()  ,tf)  =tf-t0 

s.t.  q  =  ^Clq  (3.16) 

co  =  u,  |l/|  <  1 

In  order  to  allow  a  comparison  with  published  results,  we  choose  our  initial  and 
final  conditions  as: 


X- 

=  k 

q2  q3 

q4 

co1 

co2  co3]r 

=  [o 

0  0  1 

0 

0 

or 

(3.17) 

Xf 

=  [o 

0  1  0 

0 

0 

or 

Referring  to  equations  (3.10)  and  (3.11)  these  conditions  represent  a  rest-to-rest 
rotation  maneuver  of  1 80  degrees  about  the  z-body  axis. 

As  our  first  step  in  determining  the  control  u*  that  will  drive  the  dynamic 

system  we  write  the  Hamiltonian  in  standard  form.  Since  we  have  written  the 
cost  functional  in  the  Mayer  form  it  will  not  appear  in  the  Hamiltonian.  Thus  the 
Hamiltonian  will  be  a  linear  combination  of  the  state  dynamics  and  take  the  form: 

H(X,x,u,t)  =  XT  f{x  ,u)  (3.18) 

Substituting  equations  (3.13)  and  (3.15)  into  the  Hamiltonian  equation  (3.18)  we 
have: 

X  X 

H  =  -y  (®lQ4  -  +“3^2  )  +  y  K  <h  +  -®39i  )  + 

X  X 

y(-GV72  +co2q1+co3q4)+-^(-co1q1-co2^-co3q3)+  (3.19) 

Kxu,  +  \U2  +  \>zu3 
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where  the  subscripts  on  the  Lagrange  multipliers  have  been  selected  to  aid  in 
bookkeeping. 

In  order  to  minimize  the  Hamiltonian  we  form  the  Lagrangian  of  the 
Hamiltonian  by  adjoining  the  constraint  equations  to  the  Hamiltonian  in  the  form: 

H  =  H+\iTh  (3.20) 

where  |i,are  the  KKT  multipliers  and  h(u{t),t)  is  the  control  constraint  function  in 

the  standard  form.  On  substituting  the  Hamiltonian  equation  (3.19)  and  the 
control  constraint  equation  as  defined  above  into  equation  (3.20)  the  inertial 
symmetric  problem  has  necessary  conditions, 

^,+|i  =  0  (3.21) 


with  the  complemetarity  condition, 


<0 

Ui  =  -1 

=  0 

if 

-1  <Uj<  1 

>0 

u.  =  1 

(3.22) 


The  quantity  dH/du  is  called  the  “switching  function”  in  the  literature.  The 

case  when  the  switching  function  equals  zero  for  a  non-zero  period  of  time  was 
rigorously  examined  by  Bilimoria  and  Wie  and  shown  not  to  be  time  optimal. 
Thus  we  are  left  with  a  switching  function  that  determines  when  the  optimal 
control  u*  will  switch  between  its  extreme  values.  For  this  reason  the  control 
profile  is  called  bang-bang. n 

In  order  to  validate  their  results,  Bilimoria  and  Wie  used  a  multiple¬ 
shooting  algorithm  to  solve  the  two-point  boundary  value  problem  resulting  from 
the  state  and  adjoint  equations.  This,  in  conjunction  with  the  state  and  adjoint 
equations  allowed  them  to  determine  and  evaluate  candidate  optimal  control 
solutions  and  the  resulting  trajectories.12  Here  we  employ  a  pseudospectral 
approximation  to  arrive  at  the  optimal  solution. 
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A  careful  comparison  of  previous  results  with  the  ones  obtained  here 
(Figure  2  and  Figure  3)  reveals  that  for  this  case  there  are  at  least  two  and 
potentially  four  equally  optimal  solutions.  For  the  symmetric  spacecraft  the 
precession  that  is  characteristic  of  the  optimal  solution  may  proceed  in  either 
direction  for  a  180  degree  reorientation.  A  simple  transformation  of  the  controls 
as  follows: 

Current  Previous 
ui  ->  (— 1)  x  u2 
u2  —>  L/1 
u3-^u3 

will  transform  the  current  quaternion  and  angular  rate  histories  figure  2  and 
Figure  3)  to  an  exact  match  with  those  of  Bilimoria  and  Wie.  These  transformed 
results  are  shown  in  Figure  4  and  Figure  5.  For  the  180  degree  maneuver  under 
consideration  one  could  argue  that  the  entire  maneuver  could  proceed  in  the 
opposite  direction  and  still  result  in  the  same  cost  and  hence  be  equally  optimal. 
This,  while  true,  is  precluded  by  the  quaternion  definition  convention  employed. 

Therefore,  the  results  obtained  correspond  to  the  published  results  and 
clearly  show  that  the  solution  is  not  an  eigenaxis  maneuver.  This  is  evident  from 
both  the  non-zero  quaternion  histories  of  ^  and  q2  (Figure  2)  and  the  non-zero 
angular  velocities  about  the  x  and  y  axes  (Figure  3). 
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Angular  Rates  History 


Figure  3  Inertial  Symmetric  Angular  Rate  Solution 
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Figure  5  Inertial  Symmetric  Angular  Rate  Transformed  Solution 
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The  optimal  control  solution  is  shown  in  Figure  6.  This  solution  matches 
the  previous  work  of  Bilimoria  and  Wie  and  exhibits  the  five  control  switch,  bang- 
bang  structure  that  they  observed. 


The  control  vector  obtained  is  verified  as  a  feasible  control  through  propagation 
of  the  state  dynamics.  The  initial  conditions  and  control  solution  are  used  as 
input  to  a  MATLAB®  ODE45  propagation  subroutine  which  uses  an  explicit  one- 
step  Runge-Kutta  medium  order  (4th  -  to  5th  -order)  solver^  to  verify  that  the 
control  solution  drives  the  system  from  the  given  initial  condition  to  the  desired 
final  condition.  A  linear  interpolation  was  used  to  approximate  the  control  values 
between  LGL  points.  Propagation  results  are  shown  in  Figure  7  and  Figure  8. 
The  original  solution  obtained  is  shown  in  solid  lines  overlaid  with  the  propagated 
states  shown  as  V  marks  below. 
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Figure  7  Quaternion  Propagated  Solution 


It  is  easy  to  see  that  not  only  does  the  dynamic  system  propagate  to  the  desired 
end  state  but  that  the  pseudospectral  approximation  of  the  states  closely 
matches  the  propagated  results. 


Figure  8  Angular  Rate  Propagated  Solution 
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Next  we  examine  the  necessary  conditions  for  optimality.  Recall  that 
equation  (3.21)  and  the  complementarity  conditions  of  equation  (3.22)  define  the 
switching  structure  of  the  control  vector  and  define  a  relationship  between  the 
costate  dynamics  and  KKT  multipliers.  An  inspection  of  the  switching  functions 
and  their  relationship  to  the  control  behavior  verifies  that  the  control-constraint 
pair  meet  the  KKT  conditions.  Switching  functions  for  each  axis  are  shown 
(Figure  9 ,  Figure  1 0,  and  Figure  1 1 ). 


Figure  9  Inertial  Symmetric  Spacecraft  Control  and  Switching  Function 

About  X-body  Axis 
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Figure  10  Inertial  Symmetric  Spacecraft  Control  and  Switching  Function 

About  Y-body  Axis 


Figure  1 1  Inertial  Symmetric  Spacecraft  Control  and  Switching  Function 

About  Z-body  Axis 
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Turning  to  the  Hamiltonian  equation  (3.19)  there  is  no  explicit  dependence 
on  time  and  therefore  the  Hamiltonian  will  have  a  constant  value  along  the 
optimal  trajectory. 


Recall  the  Hamiltonian  value  condition  allows  us  to  determine  the  final  value  of 
the  Hamiltonian  from  the  end-point  Lagrangian.  The  end  manifold,  in  standard 
form,  is  given  by: 


e{xf,tf) 


% 

q3f  -i 

% 

CO, 

't 

co„ 


co 


3, 


(3.23) 


Then  the  end-point  Lagrangian,  previously  defined  as: 

E  =  E+vTe  (3.24) 

allows  us  to  determine  the  final  value  of  the  Hamiltonian  as: 


H[tf]  =  _1 


(3.25) 


Thus,  the  value  of  the  Hamiltonian  will  be  -1  at  all  times  along  the  optimal 
trajectory.  The  resultant  Hamiltonian,  shown  in  Figure  12,  meets  the  necessary 
conditions  with  some  small  numerical  variation. 
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Figure  12  Inertial  Symmetric  Problem  Hamiltonian  Evolution 

The  Adjoint  equations  can  be  formed  by  differentiation  of  the  Hamiltonian. 
However,  since  the  state  variables  are  specified  at  both  the  initial  and  final 
conditions  the  adjoint  variables  will  be  free  or  unspecified  at  both  initial  and  final 
conditions.  Therefore,  the  adjoint  equations  and  terminal  transversality  of  the 
adjoint  variables  provide  no  new  information  which  will  aid  in  our  solution  to  the 
problem. 


3.  Numerical  Considerations  and  Notes 

We  have  shown  that  the  optimal  control  solution  is  feasible  and  meets  the 
necessary  conditions  derived  from  Pontryagin’s  Minimum  Principle.  Additionally, 
we  have  shown  that  the  results  obtained  closely  match  those  previously  obtained 
by  Bilimoria  and  Wie.  Next  we  compare  the  time  required  for  optimal 
reorientation  with  that  required  for  a  comparable  eigenaxis  rotation. 
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Symmetric  Maneuver  Comparison 


45 

90 

135 

180 

□  Optimal 

1.752 

2.428 

2.895 

3.254 

■  Eigenaxis 

1.777 

2.513 

3.078 

3.555 

T able  1  Comparison  of  symmetric  spacecraft  reorientation  time  for  time- 

optimal  and  eigen  axis  maneuvers 


A  comparison  of  the  time  required  for  reorientations  is  shown  in  Table  1. 
These  results  indicate  that  the  time-optimal  maneuver  is  always  faster  than  the 
optimal  eigenaxis  maneuver  except  as  noted  by  Bilimoria  and  Wie14.  The  cost 
reduction  is  shown  below  in  Table  2.  It  is  clear  that  larger  reorientation 
maneuvers  represent  a  larger  cost  benefit. 

Proper  scaling  of  numerical  problems  is  important  to  both  the  accuracy  of 
the  result  and  the  computation  time  required  in  obtaining  the  result.  No 
numerical  scaling  was  employed  in  this  algorithm.  The  nature  of  the  problem  is 
such  that  for  reasonably  large  rotational  maneuvers  ail  numerical  values  are  of 
the  same  order  of  magnitude. 
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Percent  Reduction 


Table  2  Cost  reduction  for  inertial  symmetric  time-optimal  reorientations 

over  eigen  axis  maneuvers 


4.  Conclusions 

The  published  work  by  Bilimoria  and  Wie  represents  a  landmark 
achievement  in  the  application  of  optimal  control  theory.  Their  basic  result,  that 
the  time-optimal  maneuver  is  not,  in  general,  an  eigenaxis  maneuver  will  be 
shown  to  extend  to  a  wide  variety  of  spacecraft  moments  of  inertia  and  control 
configurations. 


C.  AXISYMMETRIC  SPACECRAFT  REORIENTATIONS 

In  this  section  we  will  examine  the  time-optimal  reorientation  of 
axisymmetric  spacecraft.  Shen  and  Tsiotras15  examined  the  problem  of 
axisymmetric  reorientations  using  two  control  torques.  They  used  a  combination 
of  direct  and  indirect  methods  to  numerically  evaluate  several  representative 
maneuvers16.  We  begin  with  reorientations  about  the  axis  of  symmetry  before 

30 


examining  the  general  reorientation  case.  Finally,  we  will  examine  the 
reorientation  of  the  axis  of  symmetry  of  a  spacecraft  spinning  about  its  axis  of 
symmetry.  For  the  spinning  spacecraft  we  will  examine  the  case  where  the  spin 
rate  is  held  constant  throughout  the  maneuver.  This  amounts  to  the  two-control 
reorientation  of  Shen  and  Tsiotras.  Additionally,  we  will  examine  the  case  where 
the  spin  rate  is  allowed  to  vary  during  the  maneuver  from  a  given  initial  condition 
to  a  known  final  condition.  We  will  show  that  adding  this  third  control  torque 
results  in  the  time-optimal  maneuver  for  a  spinning  spacecraft. 

1.  Problem  Formulation 

The  inertial  axisymmetric  body  is  shown  in  Figure  13.  The  axis  of 
symmetry  is  chosen  as  the  z-body  axis  somewhat  arbitrarily  though  this  is  not  an 


Z-Body 


X-Body  Y-Body 

Figure  13  Inertial  Axisymmetric  Body 

uncommon  configuration.  As  before  the  state  of  the  spacecraft  will  be  defined 
as: 


31 


The  quaternion  kinematics  remain  as  previously  defined  in  equation  (3.12). 
However,  the  rotational  dynamics  in  the  principal  axis  frame  given  in  equation 
(3.14)  now  take  the  following  form: 

/Wi  =/xoix+(/z-/y)coycoz 

M2  =  /y°iy  +  {lx  -  4  )ff>xcoz  (3.26) 

M3  =  /zcoz 


Note  that  the  gyroscopic  term  about  the  zbody  axis  has  been  eliminated  by  the 
equal  moments  of  inertia  about  the  x  and  ybody  axes.  Thus,  the  time  rate  of 
change  of  angular  rate  is  given  by: 


co  = — 1 


M, 
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co, 


CO. 


M2 

'y 

Ms 
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(  i  -I  A 

z  y 

v  *  j 

ri -I  ' 


|coycoz 


Koxcoz 


V  7  J 


(3.27) 


These,  taken  together  with  the  quaternion  kinematics  equations  defined  earlier 
form  the  new  dynamic  constraints  on  the  axisymmetric  system.  Once  again  the 
control  parameter  is  chosen  as  external  torque  and  is  limited  to  unity. 
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2.  Time  Optimal  Maneuvers 

The  formal  statement  of  the  optimal  control  problem  is  as  follows: 


Minimize  J(x{),u{),tf)  =tf-t0 
s.t.  g  = 


“x  =  — “ 

'  \  -1  } 

z  y 

x  / 

1 

X 

V  x 

f ,  .  \ 

co 

i 

y 

\  y  j 

U3 

CO  =  — 

Z  / 

z 

|u|  <  1 

(3.28) 


a.  Non-spinning  Axisymmetric  Reorientations  About  the 
Axis  of  Symmetry 

We  will  begin  by  considering  a  family  of  prolate  spacecraft  with 
characteristics  as  given  (Table  3).  Assuming  the  spacecraft  has  uniform  mass 
distribution,  the  moments  of  inertia  about  the  principal  axes  are  given  by: 


I  y  =  I  =—mr2  +—mf2 

x  y  4  12 


(3.29) 


where,  t  is  the  length  of  the  axis  of  symmetry,  mis  the  mass  of  the  spacecraft 
and  r  is  the  radius.  The  effect  of  increasing  length  on  moment  of  inertia  is  shown 
in  Table  3. 
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Case 

Mass 

Length 

Radius 

lx 

iy 

Iz 

1 

1 

1 

1 

0.33333 

0.33333 

0.5 

II 

1 

5 

1 

2.33333 

2.33333 

0.5 

III 

1 

10 

1 

8.58333 

8.58333 

0.5 

IV 

1 

20 

1 

33.5833 

33.5833 

0.5 

V 

1 

50 

1 

208.583 

208.583 

0.5 

VII 

1 

100 

1 

833.583 

833.583 

0.5 

Table  3  Prolate  Spacecraft  Characteristics 


The  initial  and  final  conditions  will  define  a  rest-to-rest  maneuver  in  inertial  space 
representative  of  a  135  degree  rotation  about  the  axis  of  symmetry.  The  initial 
condition  will  be  defined  as  nadir  pointing  therefore  we  have  the  following: 
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Referring  to  equation  (3.18)  we  can  write  the  Hamiltonian  for  the  axisymmetric 
system  as  follows: 

X  X 

H  =  -y(©i<74  -“2^3  +w3Q2  )  +  +  «2Q4  -«3<7i  )  + 

X  X 

-y(-o \qz  +“2Qi  +“3q4)+-^(-WiQi  -©2  <*>  -©3q3)+  (3-3°) 
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z  y 
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1 

V  * 

K  x  J 

) 

l  y 

l  y  ) 

) 

where  we  have  followed  the  Lagrange  multiplier  subscript  convention  established 
earlier.  Then  minimizing  the  Hamiltonian  by  forming  the  Lagrangian  of  the 
Hamiltonian  and  evaluating  the  partial  derivative  of  the  Lagrangian  with  respect 
to  the  control  as  before  we  have, 
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\ )z  _  n 

—j—  +  (^3“0 


(3.31) 


as  necessary  conditions  for  optimal  control  solution.  Turning  our  attention  to  the 
first  prolate  spacecraft  case  (Table  3)  we  find  that  the  optimal  solution  exhibits 
characteristics  similar  to  the  symmetric  spacecraft  we  examined  previously.  The 
optimal  control  solution  is  shown  in  Figure  14.  It  displays  the  same  bang -bang 
switching  structure  that  was  previously  observed  in  the  symmetric  case. 
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Figure  14  Case  1  -  Axisymmetric  Spacecraft  Time-optimal  Control 

Solution 


The  quaternion  and  angular  rate  histories  are  shown  (Figure  15  and 
Figure  16)  and  clearly  demonstrate  that  the  time  optimal  maneuver  is  not  an 
eigenaxis  maneuver. 
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Figure  15  Case  I  -  Axisymmetric  Spacecraft  Time-Optimal  Quaternion 

History 


Figure  16  Case  I  -  Axisymmetric  Spacecraft  Time-Optimal  Angular  Rate 

History 

The  feasibility  of  the  solution  is  verified  by  propagating  the  optimal  control 
solution  through  the  state  dynamics.  The  propagation  results,  shown  in  Figure 
17  and  Figure  18,  indicate  that  the  solution  does  indeed  drive  the  state  from  the 
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known  initial  conditions  to  the  desired  final  conditions.  The  original  solution 
obtained  is  shown  in  solid  lines  overlaid  with  the  propagated  states  shown  as  V 
marks  below. 


Figure  17  Axisymmetric  Spacecraft  Quaternion  Solution  Validation  by 

Propagation 


Figure  18  Axisymmetric  Spacecraft  Angular  Rate  Solution  Validation 
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The  necessary  conditions  for  optimality  that  we  can  evaluate 
directly  include  the  switching  structure  obtained  from  the  Hamiltonian 
minimization  (Equation  (3.31)),  and  the  behavior  of  the  Hamiltonian  over  time. 
The  switching  structure  is  shown  in  Figure  19  through  Figure  21.  The  control 
solution  has  been  overlaid  to  further  illustrate  the  relationship  between  the 
switching  function  and  the  control  solution. 
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Figure  19  Axisymmetric  Spacecraft  x-axis  Switching  Structure  and 

Control 
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Figure  20  Axisymmetric  Spacecraft  y-axis  Switching  Structure  and 

Control 


Figure  21  Axisymmetric  Spacecraft  z-axis  Switching  Structure  and 

Control 
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Inspecting  the  Hamiltonian  for  the  axisymmetric  case  (Equation 
(3.30))  reveals  no  direct  dependence  on  time.  Therefore  we  can  see  that  the 
Hamiltonian  should  be  a  constant  over  the  interval  under  consideration. 
Additionally,  forming  the  end -state  Lagrangian  as  before  in  equation  (3.24)  gives 
the  final  value  of  the  Hamiltonian.  The  evolution  of  the  Hamiltonian  over  time  is 
shown  in  Figure  22. 


Figure  22  Case  1  -  Axisymmetric  Hamiltonian  Evolution  and 

Transversality 

Our  analysis  of  the  solution  indicates  that  it  is  a  feasible  solution  to 
the  time-optimal  reorientation  problem.  Additionally,  the  solution  meets  the 
necessary  conditions  for  optimality  derived  from  Pontryagin’s  Minimum  Principle. 

Now  examine  the  effects  of  increasing  the  length  of  the  symmetry 
axis  of  a  constant  mass  prolate  spacecraft.  From  equations  (3.29)  and  Table  3 
we  can  see  that  the  moment  of  inertia  about  the  axis  cf  symmetry  remains 
constant  while  the  remaining  moments  of  inertia  grow  proportionally  to  length 
squared.  This  is  shown  in  Figure  23.  These  increasing  moments  of  inertia  have 
two  effects.  First,  the  angular  acceleration  about  the  x  &  y-body  axes  is  reduced 
proportional  to  the  increase  in  moment  of  inertia.  Physically,  this  is  a  decrease  in 
the  ability  of  a  constant  torque  to  cause  an  angular  acceleration  of  an  increasing 
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mass  moment  of  inertia.  However,  in  addition  to  the  decreasing  control 
effectiveness,  we  see  that  the  switching  functions  of  the  controls  about  the  x  &  y- 
axes  are  numerically  vanishing  due  to  the  increasing  moment  of  inertia.  The  y- 
axis  switching  function  for  case  III  illustrates  this  effect  and  is  shown  in  Figure  25. 
Increasing  the  value  of  the  control  torque  available  is  not  sufficient  to  counter  the 
vanishing  switching  function. 


Moment  of  Inertia 


Length 


Figure  23  Constant  Mass  Prolate  Spacecraft  Moment  of  Inertia  versus 

Length  of  Symmetry  Axis 

Therefore  the  time-optimal  reorientation  maneuver  about  the  symmetry  axis 
approaches  an  eigenaxis  maneuver  as  the  length  of  the  body  grows  in  a  constant 
mass  spacecraft. 
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Figure  24  Axisymmetric  Reorientation  Maneuvers  About  Symmetry  Axis 

versus  Symmetry  Axis  Length 

It  is  interesting  to  note  that  the  structure  of  switching  function  is  not 
significantly  affected  by  the  increasing  moments  of  inertia.  The  switching 
function  for  the  prolate  spacecraft  of  case  III  (Table  3),  where  the  length  of  the 
symmetry  axis  is  10  times  longer  than  the  original,  displays  the  same  switching 
structure  as  the  original.  The  magnitude  is  however,  reduced  by  several  orders 
of  magnitude.  The  switching  function  for  the  y-axis  control  of  Case  III  is  included 
(Figure  25)  for  comparison  to  that  of  Case  I  (Figure  20). 
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Figure  25  Prolate  Axisymmetric  Spacecraft  Case  III  -  Switching  function 

b.  Non-spinning  Axisymmetric  Reorientations  of  the  Axis 
of  Symmetry 

The  case  of  non-spinning  axisymmetric  spacecraft  reorientations  of 
the  axis  of  symmetry  is  mathematically  no  different  from  the  previous  section. 
The  problem  formulation  and  necessary  conditions  for  optimality  remain 
unchanged.  The  optimal  maneuvers  display  the  same  characteristic  precession 
about  the  eigenaxis  and  involve  control  torques  about  all  three  axes.  It  is 
interesting  to  note  that  the  control  switching  structure  behaves  in  a  manner  very 
similar  to  that  observed  by  Bilimoria  and  Wie  in  their  work  on  symmetric 
reorientations.17  That  is,  for  small  angle  reorientations  we  observe  a  sequential 
seven-switch  structure  and  for  large  angle  maneuvers  we  observe  a  sequential 
five-switch  structure.  The  control  solutions  for  two  representative  maneuvers  are 
shown.  In  Figure  26,  we  see  the  control  solution  for  a  large  angle  maneuver 
about  the  x-axis.  The  representative  maneuver  is  chosen  as  a  135  degree 
rotation.  The  five-switch-sequential  structure  is  clearly  evident.  In  Figure  27,  the 
control  solution  for  a  representative  small  angle  maneuver  about  the  x-axis  is 
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Figure  26  Axisymmetric  Spacecraft  Optimal  Control  Solution  for  X-axis 
Large  Angle  Rotation  (135  Degree  Rotation) 


Figure  27  Axisymmetric  Spacecraft  Optimal  Control  Solution  for  X-axis 
Small  Angle  Rotation  (60  Degree  Rotation) 
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shown.  The  maneuver,  a  60  degree  rotation  results  in  a  seven-switch-sequential 
solution.  This  suggests  that  there  is  a  dividing  point,  as  Bilimoria  and  Wie 
observed  for  the  symmetric  case,  where  the  time-optimal  solution  changes  from 
the  seven-switch  to  the  five-switch  structure. 


c.  Reorientations  of  the  Spinning  Axis  of  Symmetry 

Next  we  consider  the  reorientation  of  a  spinning  axisymmetric 
spacecraft.  Shen  and  Tsiotras  also  considered  this  case  where  the  rigid  body 
was  subject  to  only  two  control  torques  which  spanned  the  plane  perpendicular  to 
the  axis  of  symmetry.  They  concluded  that  two  torques  were  sufficient  to  achieve 
a  time-optimal  maneuver.18  In  this  section  we  will  show  that  a  third  torque  about 
the  symmetry  axis,  if  available,  further  reduces  the  objective  function  and  is  the 
true  time-optimal  solution. 

Consider  the  axisymmetric  spacecraft  of  Figure  13,  which  is 
spinning  about  the  z-axis,  the  axis  of  symmetry.  Euler’  equations  remain 
unchanged  from  equation  (3.26)  and  are  repeated  here  for  the  convenience  of 
the  reader: 

=/xojx+(/z-/y)cor(oz 
M2  =/rojy  +  (/x-/z)coxcoz 
M3  =  ip 6Z 

Using  the  formulation  of  Shen  and  Tsiotras  suggested  by  Tsiotras  and  Longuski19 
the  orientation  of  the  h3  inertial  axis  of  the  inertial  frame  given  by 

n  =  (hv n2, 4)  with  respect  to  the  body  frame  can  be  represented  by  two  variables 
defined  as: 

P  -a 

w \  =  — —  w„  = - 

1+Y  1+y 

where  w1  and  w2  obey  the  differential  equations: 
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•  UJ  /  p  2  \ 

=  CO zw2+  CO ywj/z2  +  —  (1  +  i/v,  -W2J 


24 

2 

CO, 


w2=-wzw,+0)xw,w2+^(l 


+  w\-  wf 


Readers  are  directed  to  the  references  [18,19]  for  a  complete  derivation  of  this 
parameterization.  Using  this  formulation  we  can  state  the  time-optimal 
reorientation  of  the  spin  axis  as  follows: 


Minimize  J(x{),u{),tf)  =  tf-t0 
s.t. 
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+  coy  w,w2  +^(1  +  Wf- w22 ) 

CO 


w2  =  -CO  Zw1  +  (OxW,W2  +— (l 


+  w2  -  wf 


) 


lulsl 


Following  the  numerical  example  of  Shen  and  Tsiotras  we  establish  the  following 
spacecraft  parameters: 


4  =  4 
4=4 

4=2 


The  Hamiltonian  and  necessary  conditions  are  available  in  the  reference  and  not 
repeated  here.  Using  initial  and  final  conditions  given  as: 
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x  =  [co1  co2  co3  w1  w2f 
x0  =  [0  0  -0.5  1.5  -0.5]r 
xf  =  [0  0  -0.5  0  Of 

we  can  duplicate  the  results  of  Shen  and  Tsiotras  by  setting  the  torque  M3  =  0. 
The  results  are  shown  in  Figure  28  and  Figure  29.  The  maneuver,  defined  as  a 
115.38  degree  reorientation  of  b3  to  h3,  and  is  completed  in  2.6142  seconds. 
Published  results  indicated  a  minimum  maneuver  time  of  2.61  seconds.20 


Figure  28  Angular  Rate  History  for  Constant  Spin  Rate  Time-optimal 

Maneuver 
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Figure  29  Time  Histories  of  Wi  and  W2  for  Constant  Spin  Time-optimal 

Maneuver 


However,  if  the  control  u3  is  available,  and  the  boundary  conditions 

are  enforced  such  that  the  spin  rate  is  allowed  to  vary  throughout  the  maneuver, 
then  the  true  time-optimal  solution  is  found. 


Figure  30  Angular  Rate  History  for  Time-optimal  Spinning  Reorientation 

Maneuver 
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If  spin  rate  of  the  z-axis  is  allowed  to  vary  the  time-optimal 
maneuver  solution  contains  a  significant  change  in  spin  rate  as  shown  in  Figure 
30.  The  improved  time  for  maneuver  completion  is  now  2.513  seconds.  This  is  a 
reduction  of  3.87%  from  the  previously  assumed  time-optimal  solution.  The 
three-axis  control  time-optimal  w solution  is  shown  in  Figure  31 . 

The  control  solution  obtained  (Figure  32)  is  bang-bang  in  all  three 
axes.  As  before  the  control  solution  is  propagated  through  the  state  dynamics  to 
verify  that  the  solution  is  feasible.  The  original  solution  obtained  is  shown  in  solid 
lines  overlaid  with  the  propagated  states  shown  as  V  marks  below.  It  is  clear 
(Figure  33  and  Figure  34)  that  the  spacecraft  states  properly  propagate  from  the 
given  initial  conditions  to  the  desired  final  conditions. 


Figure  31  Time  Histories  of  w1  and  W2  for  Spinning  Spacecraft  Time- 

optimal  Maneuver 
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Figure  32  Spinning  Axisymmetric  Spacecraft  Time-optimal  Control 

Solution 


Figure  33  Spinning  Axisymmetric  Spacecraft  Angular  Rate  Solution 

Validation  by  Propagation 
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Figure  34  Spinning  Axisymmetric  Spacecraft  W  History  Solution 

Validation  by  Propagation 

Minimization  of  the  Hamiltonian  with  respect  to  the  control  vector 
allows  us  to  establish  the  switching  functions.  These  are  shown,  overlaid  with 
the  unity  scaled  control  solution  (Figure  35,  Figure  36,  and  Figure  37). 
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Figure  35  Spinning  Axisymmetric  Spacecraft  X-axis  Switching  Function 

and  Normalized  Control 

Finally,  we  evaluate  the  Hamiltonian  and  observe  it  is  constant  over 
time  and  numerically  equal  to  -1.  This  is  expected  for  the  Mayer  formulation  of 
the  cost  function  when  the  Hamiltonian  has  no  direct  time  dependence. 
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Figure  37  Spinning  Axisymmetric  Spacecraft  Z-axis  Switching  Function 

and  Normalized  Control 


Figure  38  Spinning  Axisymmetric  Spacecraft  Time-optimal  Solution 

Hamiltonian 
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3.  Numerical  Considerations  and  Notes 

Operationally,  it  is  often  preferable  to  maintain  a  constant  satellite  spin 
rate  about  the  axis  of  symmetry.  However,  we  have  shown  that  if  a  third  torque 
is  available  about  the  axis  of  symmetry  then  the  time-optimal  solution  will  involve 
a  change  in  the  spacecraft  spin  rate.  This  change  in  spin  rate  is  appreciable  and 
results  in  a  measurable  time  savings  over  the  two  control  torque  solution. 

4.  Conclusions 

In  this  section  we  examined  the  time-optimal  maneuver  for  an 
axisymmetric  spacecraft  with  three  independent  control  torques.  We  saw  that  for 
maneuvers  about  the  axis  of  symmetry  the  relative  moment  of  inertias  have 
significant  effects  on  the  maneuvers.  As  the  moment  of  inertia  about  the 
symmetry  axis  increased  the  maneuver  approached  the  eigenaxis  maneuver  in 
both  time  required  and  spacecraft  response.  Additionally,  we  demonstrated  that 
while  reorientation  of  the  spin  axis  is  possible  with  two  control  torques  spanning 
the  plane  perpendicular  to  the  spin  axis,  the  addition  of  a  third  control  torque 
about  the  axis  of  symmetry  further  reduces  the  objective  function  and  is  the  true 
time-optimal  solution.  This  third  control  torque  is  in  general  assumed  to  be 
available  as  it  was  required  to  generate  the  spinning  motion. 

D.  ASYMMETRIC  SPACECRAFT  REORIENTATION  MANEUVERS 

In  this  section  we  will  numerically  investigate  the  time-optimal  reorientation 
of  a  rigid  asymmetric  body.  Livenh  and  Wie21  presented  an  extensive  analytical 
analysis  of  the  asymmetric  reorientation  problem  under  constant  body-fixed 
torques.  Additionally,  the  work  of  Proulx  and  Ross22  determined  an  admissible 
switching  structure  which  was  elegantly  illustrated  by  the  traversal  of  a  unit  cube. 
Using  this  to  limit  the  search  space  a  combination  of  a  genetic  algorithm  and 
pseudospectral  method  was  used  to  obtain  the  optimal  solution.  Additionally, 
they  suggested  a  method  of  evaluating  the  “optimality”  of  a  solution  by  evaluating 
the  Hamiltonian  derived  from  the  costates  obtained  through  the  Covector 
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Mapping  Theorem.  This  method  of  evaluating  compliance  with  the  Minimum 
Principle  is  employed  throughout  this  work. 


1.  Problem  Formulation 

In  order  to  facilitate  future  work  and  extension  to  other  control  actuators 
the  orientation  of  the  asymmetric  spacecraft  will  be  represented  with  respect  to 
an  Earth  centered  inertial  reference  frame.  This  will  require  the  incorporation  of 
orbital  velocity  and  reference  frame  transformations  to  properly  represent  the 
dynamic  constraints. 

As  before,  for  a  rigid  body  the  applied  torque  about  the  center  of  mass  is 
equal  to  the  time  rate  of  change  of  the  angular  momentum.23 


From  the  transportation  theorem  we  have, 


d 

dt 


+ 


VxH 


Following  the  previous  development  we  obtain  Euler’s  equations: 

=/x«x+(4-/yK0)z 

M2  =  ljby  +  (lx  -  4  K®z 

M3  =  /zc6z  +  (ly  -lx) G)xG)y 


(3.32) 


The  quaternion  kinematics  equations  are  unchanged  from  previous  (equation 
(3.13))  except  that  emphasis  is  placed  on  the  notation. 


4  =^[coi<74  -®2<73  +0)3<72] 

4  =^[«lQ3+^Q4-«3Ql] 
4  =^[-°V72  +“2Qi  +“3^] 
Q4=li-(»a-(»zQ2-^ch] 


(3.33) 
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As  shown  above  the  angular  rates  (co1,co2,co3 )  are  of  the  body  with  respect 

to  the  orbit  frame.  For  convenience  we  choose  to  define  the  state  of  the 
spacecraft  in  terms  of  the  quaternion  vector  and  the  angular  rates  with  respect  to 
the  Newtonian  frame.  Thus  the  state  vector  is  given  by: 


2.  Time  Optimal  Maneuvers 

The  formal  statement  of  the  optimal  control  problem  is  as  follows: 

Minimize  J(x{),u{),tf)  =tf-t0 

s.t.  q  =  —  £lq 
-  2  - 


u,  <1 


The  maneuver  under  consideration  will  be  rest-to-rest  in  the  orbit  frame; 
however,  the  final  angular  velocity  in  the  inertial  frame  will  depend  on  the  final 
attitude.  This  is  clear  from  equation  (3.8).  The  maneuver  under  consideration  is 
an  x-axis  rotation  typically  135  degrees.  This  maneuver  magnitude  was  selected 
based  on  the  anticipated  orbital  parameters  of  NPAST1,  the  potential  test  bed  for 
later  algorithms.  The  initial  condition  was  chosen  somewhat  arbitrarily  as  nadir 
pointing.  Then  the  problem  initial  and  final  conditions  may  be  presented  in 
standard  form  as  follows: 
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x  =  [q i  q2  q3  qA  cox  coy  coj 

Xo=[0  0  0  1  0  -co0  0]r 


£  = 

sin 

( 1 35°  ^ 
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0 

0  cos 

f  1 35°  ^ 

p 

G>xU)  “y(0  ®zU) 

- 

v  J 

v  J 

- 

Additionally,  spacecraft  moment  of  inertias  for  subsequent  examples  have  been 
selected  to  match  the  planned  NPSAT1  moment  of  inertias  and  are  given  as 
follows: 

lx=  5  kg.m2 
ly  =  5.1  kg*m2 
/z  =  2  kg*m2 

As  before,  the  first  step  in  our  approach  is  to  form  and  minimize  the 
Hamiltonian.  Recall  that  the  Hamiltonian  is  a  function  of  the  state,  control  and 
Lagrange  multipliers. 

H(X,x,u,t)  =  XTf{x,u) 

Since  quaternion  kinematics  are  typically  written  in  terms  of  the  body  angular 
rates  with  respect  to  the  orbit  frame,  a  lengthy  algebraic  process  of  coordinate 
transformations  is  required  to  properly  write  the  Hamiltonian.  Once  completed 
the  Hamiltonian  takes  the  following  form: 


X  X 

H  =  (°>x<74  -°V? 3  +“zQ2  +“oQ3  )  +  yK<73  +CV74  -G>z<7l  +“0^4  ) 

X  X 

+  -y ( -™xq2  +“yQi  +co2Qz  -<o0Qi )  +  -^(-°V7i  - “yQ2 - «zq3 - “0Q2) 


(3.34) 
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In  accordance  with  the  Minimum  Principle  the  Lagrangian  of  the  Hamiltonian  is 
formed  and  partial  derivatives  with  respect  to  the  control  vector  are  formed  which 
establish  the  control  switching  functions. 


an 

du^ 

})H 

du2 

})H 

du3 
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L 

K 

'  / 


-+  Pi  -  0 


■  +  p2  —  o 


■—j—  +  M-3  -0 


(3.35) 


The  reader  should  note  that  the  switching  functions  of  the  asymmetric  case  are 
no  different  than  those  of  the  axisymmetric  case  given  in  equation(3.31).  This 
would  correctly  lead  us  to  surmise  that  as  the  moment  of  inertias  approach  those 
previously  studied  the  results  would  closely  match  those  previously  obtained. 

The  control  solution  for  the  asymmetric  body  under  consideration  is  shown 
in  Figure  39.  It  clearly  exhibits  the  structure  we  have  come  to  expect.  In  general, 
the  control  is  bang-bang  in  all  three  axes.  There  are  five  switches  in  all, 
characteristic  of  a  large  angle  slew  with  a  single-switch  in  the  primary  axis  of 
rotation. 
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Asymmetric  Control  History 
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Figure  39  Asymmetric  Spacecraft  Time-optimal  Maneuver  Control 

Solution 

The  process  of  solution  validation  begins  with  propagating  the  candidate 
solution  through  the  state  dynamics.  A  feasible  solution  must  drive  the 
spacecraft  from  its  known  initial  state  to  the  desired  end  state.  The  calculated 
state  histories  are  shown  figure  40  &  Figure  41)  with  the  propagation  results 
(Figure  42  &  Figure  43).  The  original  solution  obtained  is  shown  in  solid  lines 
overlaid  with  the  propagated  states  shown  as  V  marks  below. 
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Figure  40  Asymmetric  Spacecraft  Time-optimal  Maneuver  Quaternion 

History 


Figure  4 1  Asymmetric  Spacecraft  Time-optimal  Maneuver  Angular  Rate 

History 
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Figure  42  Asymmetric  Spacecraft  Quaternion  Solution  Validation  by 

Propagation. 


Figure  43  Asymmetric  Spacecraft  Angular  Rate  Solution  Validation  by 

Propagation. 
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The  solution,  confirmed  feasible,  is  evaluated  for  optimality  by  observing 
the  switching  functions  and  Hamiltonian  evolution.  The  switching  functions, 
equations  (3.35)  are  shown  graphically  with  the  overlaid  control  solution. 


Figure  44  Asymmetric  Spacecraft  X-axis  Switching  Function  and  Control 

Solution 


Figure  45  Asymmetric  Spacecraft  Y-axis  Switching  Function  and  Control 

Solution 
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Figure  46  Asymmetric  Spacecraft  Z-axis  Switching  Function  and  Control 

Solution. 

These  figures  illustrate  both  the  Hamiltonian  minimization  and  the  system’s 
compliance  with  the  KKT  conditions. 

Inspection  of  the  Hamiltonian  reveals  no  direct  dependence  on  time.  Thus 
we  expect  a  constant  value  Hamiltonian.  The  final  value  of  the  Hamiltonian  is 
given  by  the  transversality  condition,  equation(3.25).  Taken  together  we  expect 
a  constant  Hamiltonian  of  -1.  Figure  47  shows  the  Hamiltonian  with  the 
optimality  characteristics  predicted. 

Finally,  we  compare  the  time-optimal  solution  with  the  eigenaxis  maneuver 
once  theorized  as  nearly  time-optimal.  The  results,  (Table  4)  show  a  significant 
reduction  in  maneuver  time. 
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Asymmetric  Maneuver  Comparison 


Eigenaxis 

6.8742  seconds 

Time-optimal 

6.1737  seconds 

Reduction 

10.19% 

T able  4  Asymmetric  Maneuver  Comparison 


Figure  47  Asymmetric  Spacecraft  Time-optimal  Hamiltonian  Evolution 

and  Transversality 

3.  Numerical  Considerations  and  Notes 

The  problem  was  transformed  into  the  orbital  frame  as  a  stepping  stone  to 
the  magnetic  torque  problem.  Orbital  position  and  velocity  will  be  necessary  for 
magnetic  field  computations.  Additionally,  this  problem  formulation  incorporated 
linear  scaling  factors.  While  not  strictly  necessary  for  the  problem  formulation 
under  consideration,  incorporating  scaling  factors  in  the  independent  torque 
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problem  allowed  validation  of  the  scaling  algorithm  prior  to  the  incorporation  of 
the  more  complicated  magnetic  field  calculations. 

For  evaluation,  the  torque  available  was  reduced  form  1  Newton-meter  to 
0.01  Newton-meters.  Without  scaling,  it  was  noted  that  costate  estimates 
increased  by  4  orders  of  magnitude.  In  this  work  we  have  found  that  large 
costate  estimates  have  generally  been  an  indication  of  a  numerically  poor 
problem  formulation.  Results  in  these  cases  have  generally  displayed  poor 
accuracy  and  long  computational  times.  Numerical  scaling  was  implemented  in 
the  following  form: 

u  =  kmu 
co  =  kjD 

T  =  ktt 

c/co  _  c/co  (  km  ' 

dt  dt  k, 

v  1  J 

where  the  overbar  indicates  a  scaled  variable. 

4.  Conclusions 

In  this  section  we  have  shown  that  the  asymmetric  problem  displays  many 
of  the  characteristics  previously  noted  in  symmetric  and  axisymmetric 
configurations.  The  time-optimal  solution,  as  we  might  have  predicted,  is  not  an 
eigenaxis  maneuver  but  instead  is  bang-bang  in  all  control  axes.  The  sequential 
switching  structure  theorized  and  observed  by  Proulx  and  Ross24  was  observed 
for  the  problem.  The  number  of  possible  configurations  for  asymmetric 
spacecraft  combined  with  the  possible  slew  maneuvers  limits  the  extent  to  which 
numerical  analysis  can  be  used  to  form  general  conclusions.  However,  the 
method  employed  allows  us  to  generate  the  time-optimal  solution  to  the 
asymmetric  configuration,  validate  the  feasibility  of  the  candidate  solution,  and 
evaluate  its  compliance  with  the  Minimum  Principle. 
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IV.  MAGNETIC  TORQUE  CONT  ROL 


A.  INTRODUCTION 

Magnetic  torque  control  has  been  effectively  used  for  momentum 
management  in  zero-momentum  systems  on  many  spacecraft  systems.  The 
methods  are  well  known  and  flight  tested.  The  use  of  magnetic  torque  control  for 
spacecraft  three-axis  stabilization  is  less  well  known  but  the  body  of  research  is 
growing.  Magnetic  torque  control  represents  a  low  cost  method  to  control  small 
spacecraft  in  reasonably  low  earth  orbit.  In  this  section  we  examine  the  basics  of 
magnetic  torque  control  and  the  solution  to  the  time-optimal  slew  of  spacecraft 
using  magnetic  torque  generating  devices.  This  problem  is  significantly  more 
complicated  then  the  idealized  actuator  problem  of  Chapter  III.  This  is  due  to  the 
resultant  cross-product  torque  generation  and  varying  magnetic  field.  It  is  well 
known  that  there  are  body-frame  orientations  where  no  torque  can  be  generated 
in  specific  directions. 

Junkins  and  Turner,  reference  [1],  discuss  the  magnetic  time-optimal 
control  of  spin-stabilized  spacecraft.  They  were  able  to  solve  the  open  loop 
problem  for  spin-axis  reorientation  and  implement  their  solution  on  the  NOVA-1 
spacecraft  in  198 1.1 

B.  BASIC  MAGNETIC  TORQUE  ATTITUDE  CONTROL 

The  magnetic  moment  generated  within  the  spacecraft,  whether 
generated  intentionally  or  inadvertently,  interacts  with  the  Earth’s  magnetic  field 
to  produce  a  torque  according  to: 

fB  =  mxB 

where  m  is  the  magnetic  dipole  moment  generated  inside  the  spacecraft  body 
and  B  is  the  Earth’s  magnetic  field  intensity.2 

*  Magnetic  dipole  moment  is  represented  as  a  lower  case  m  in  order  to  distinguish  it  from  an 
applied  torque  which  is  commonly  represented  as  an  upper  case  M. 
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The  Earth’s  magnetic  field  intensity  as  approximated  by  McElvain  (1962) 


is: 


Bx(t) 

By(t) 

Bz(t) 


cos(co0f)sin(/m) 
-cos  (/J 
2sin(co0t)sin(/m) 


(4.1) 


where,  im  is  the  inclination  of  the  satellite  orbit  with  respect  to  the  magnetic 
equator,  R  is  the  semi -major  axis  of  the  orbit,  co0  is  the  orbital  angular  velocity,  t 
is  zero  at  the  point  in  the  satellite  orbit  where  the  ascending  node  crosses  the 
equator,  and  is  the  magnetic  dipole  strength  of  the  Earth  (circa  1 975)  given  as: 

\l,  =  7.96x1 0'5Wb-m 

So  we  see  that  magnetic  field  intensity  decreases  rapidly  with  orbital 
altitude.  A  typical  magnetic  field  approximation  is  shown  for  a  satellite  in  a 
circular  orbit  with  inclination  of  35.4  degrees  at  altitude  of  560  Km.  As  expected 
the  behavior  in  the  x-z  planes  is  harmonic  (Figure  48)  at  the  orbital  frequency. 

The  magnetic  model  used  for  this  problem  differs  from  equation  (4.1).  The 
above  model  is  reasonably  accurate  as  a  first-order  approximation  but  does  not 
take  into  account  the  rotation  of  the  earth.  For  this  reason,  a  model  with  slightly 
higher  fidelity  was  adopted.  Shown  below,  as  equation  (4.2),  is  a  model  adopted 
from  Wertz,  Spacecraft  Attitude  Determination  and  Control.3  This  model 
assumes  no  orbit  precession  but  does  allow  for  the  Earth’s  rotation. 
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A" 

- 

K 

B 

Vo 

BZo 

- 

cos(a)  [cos(e)sin(/)-sin(e)cos(|cos(L/)]-sin(a)sin(e)sin(L/) 
-cos(e  )cos  (/)  -  sin(e  )sin(/)sin(L/) 
2sin(a)[cos(£)sin(/')-sin(£)cos(i)cos(tv)]  +  2cos(a)sin(£)sin(t7) 


where, 


K  =  7.943x1015-^- 
amp 

R  =  radius  from  Earth  c.m.  to  s/c  c.m.  (m) 
a=co  ot 

u  =  K>et ,  coe  =  Earth  rotational  frequency  (rad/sec) 
£  =  1 1 .7°  magnetic  dipole  tilt 
/=  orbit  inclination 


Figure  48  Earth  Magnetic  Field 


(4.2) 
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Recall  Euler’s  equations  (repeated  from  previous  for  the  convenience  of 
the  reader). 

=/x«x+(/z-/y)°)y0)z 

M2  =/ycoy  +  (/x-/z)o)xo)z 
M3  =/zcoz  +  (/y 

The  external  torque  can  now  be  replaced  by  the  torque  generated  by  the 
interaction  of  the  Earth’s  magnetic  field  and  the  satellite  magnetic  dipole  moment. 
The  rotational  dynamics  equation  of  motion  in  vector  form  then  becomes, 

/(6+  cox/co=  mxB  (4.3) 

where  angular  velocities  are  by  necessity  referenced  to  the  Newtonian  frame 
and, 

m2Bz-m:,By 

mxB  =  m3Sx-m1Sz  (4.4) 

m,B  -  m^B 

As  written,  the  cross-product  of  equation  (4.4)  is  meaningless  since  the 
dipole  moments  (m,)  are  in  the  body  frame  and  the  magnetic  field  components 

are  represented  in  the  orbit  frame.  Therefore,  the  components  of  the  magnetic 
field  which  are  in  the  orbit  frame  must  be  rotated  into  the  spacecraft  frame  by: 

Bb  =  BC°B0  (4.5) 

where  the  DCM  was  previously  defined  in  equation  (3.5).  Then  the  magnetic 
field  in  the  body  frame  in  expanded  form  is  given  by: 

ei  =  Bx0 (  4  -  4  -  4  +  Ql )  +  2By0  ( 1 q^2  +  Q3Q4 )  +  2BZo  (q,q3  +  q2 q 4 ) 

S2  =  2 BXq {qq2-q3q4)  +  Byo{4- cfi -c?3+q24)  +  2BZo{ q2q3  +  q1q4 )  (4.6) 

B3  =  2 BXo(q,q3  +  q2q4 )  +  2 BJq2q3  -  q,q4 )  +  BZq (q23  -  q2  -  q\  +  q\) 
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Finally,  by  substituting  equation  (4.6)  into  equation  (4.3)  we  can  form  the 
dynamic  constraint  equations  as  follows, 


“x=f(m2e3-m3e2)  +  /C1 

X 

«v  =  j{mA  -  m,a, )  +  /c2coxcoy, 

y 

“z=f(mie2-m2ei)  +  ^3 


k  ly~  4 

X 

k  =LzLl 

2  ly 

k  =k± 

3  /, 


These  represent  the  rotational  dynamical  equations  of  motion  for  our  spacecraft, 
Earth  magnetic  field  system. 


C.  TIME-OPTIMAL  MAGNETIC  TORQUE  CONTROLLED  SLEW 

In  this  section  we  consider  the  time-optimal  reorientation  of  a  spacecraft 
with  magnetic  torque  control.  The  maneuver  is  defined  as  rest-to-rest  in  the  orbit 
frame  where  the  initial  and  final  states  are  given. 


1.  Problem  Formulation 

The  spacecraft  state  is  defined  by  its  position  and  angular  velocity. 


x  ■■ 


q 

to 


The  position  is  represented  by  a  four-element  quaternion  vector,  we  have 
previously  adopted  the  convention  that  the  fourth  element  of  the  quaternion 
vector  is  the  scalar  quantity.  The  angular  velocity  is  in  body  coordinates  with 
respect  to  the  Newtonian  frame. 

We  choose  our  control  parameter  as  the  magnetic  dipole  moment  of  the 
torque  rods.  Dipole  moment  is  controlled  by  current  flow  however  the  response 
of  dipole  moment  to  changes  in  current  can  be  considered  instantaneous.  Since 
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the  actual  torque  generated  is  limited  by  the  maximum  dipole  moment  of  the 
torque  rod,  we  impose  a  bound  on  the  control  dipole  moment. 


u  = 


m2 

m3 


->  |n^|<  30 A  m2 


(4.7) 


The  minimum  time,  rest-to-rest,  reorientation  problem  may  then  be  stated  as 
follows: 

Determine  the  controls  [tv/,  u2*,u3  *] that  drive  the  spacecraft  from  its  initial 
rest  position,  given  by  [Xq] to  its  final  rest  position  given  by  [xf] while  minimizing 
the  cost  function: 

J(x(-),u(-),tf)  =  tf-t0 

where  we  have  used  the  Mayer  form  of  the  cost  function,  subject  to  the  following 
constraints: 

Control  Constraint:  The  control  constraint  is  defined  in  the  standard  form, 

hL{t)<h{u,t)<hu{t) 

then  our  control  constraint  can  be  written, 

-30  <  u,  <  30 
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Dynamic  Constraints: 


^1=^KC?4-(02C?3+(03^] 

4  =^[“lQ3+^Q4-«3Ql] 
4  =^[-°)1C?2  +“2Qi  +“3^4] 
4  =^[-®iQi-^^-«3q3] 


m2B3  ~  m3e2)  + 

X 

_ N 

1 

_ 

II 

4sT 

“y  =  y-(m3S1  -  mA  )  + 

y 

/  -  / 

>< 

N 

il 

CM 

“z  =  y-( mA  -  m2B, )  +  /c3coxcoy , 

/  -  / 

k=  y 

3  L 

(4.8) 


2.  Solving  the  Optimal  Control  Problem 

The  first  step  in  solving  the  optimal  control  problem  is  to  form  the 
Hamiltonian.  Basic  format  of  the  Hamiltonian  is  repeated  heret, 

H(X,  x,  u,  t)  =  F(x,  u,  t)  +  XTf(x,  u ;  t) 

Since  the  cost  functional  was  formulated  without  a  Lagrange  cost  term  the 
Hamiltonian  reduces  to  the  following. 

H{rk,x,u,t)  =  XTf(x,u,t ) 

Again,  since  the  quaternion  kinematics  are  written  in  terms  of  angular  velocity  of 
the  body  with  respect  to  the  orbit  frame,  equation  (4.8),  a  lengthy  algebraic 
rotation  sequence  is  required  to  write  the  Hamiltonian  in  standard  form: 


t  Recall  that  in  this  notation  F  is  the  Lagrange  (running)  cost  and  f  are  the  state  dynamics. 
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H  =  Y^xq^  _co^  +co^2  +C0oQ3)  +  -|lKQ3  +“yq4-wzQ1  +  co0Q4) 

X  X 

+  -^(-«xQ2  +®yQl  +“zQ2-W0Ql)  +  -^(-®xQl  -<V72  -«zC?3-W0Q2)  + 

-  U3B2)  +  M>y<°z}  +  V  |f  (U3ei  -  Uie3  )  + 

^coz|y-(Wie2-^2^)  +  /f3«x®y 

The  subscripts  on  the  Lagrange  multipliers  have  been  chosen  for  bookkeeping 
purposes.  The  control  vector  is  defined  in  equation  (4.7). 

Now,  according  to  Pontryagin’s  principle  the  control  which  minimizes  the 
cost  functional  must  meet  the  conditions  we  established  earlier.  It  is  however, 
important  to  note  that  not  all  of  these  conditions  reveal  usable  information  about 
the  nature  of  the  problem. 


a.  Hamiltonian  Minimization 

We  know  from  previous  work  that  a  necessary  condition  for  the 
Hamiltonian  to  be  a  minimum  with  respect  to  the  control  variable  is  that  the 
partial  derivative  of  the  Lagrangian  with  respect  to  the  control  equals  zero.  In  this 
case  the  control  that  satisfies  this  condition  must  also  lie  within  the  control 
constraint  set.  We  apply  Lagrange  multipliers  in  the  form, 

H(X,x,u,t,[i)  =  H(k, x,  u,  t)  +  \xTh{u ,  t)  (4.10) 

Then,  by  inspection  of  equations  (4.9)  and  (4.10),  we  can  write  the 
necessary  conditions  for  the  Hamiltonian  minimization  as  follows. 
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+  m  =o 

4  2  4 

^B3-^B1+p2=0  (4.11) 

X  z 

^^-^62  +  1X3  =0 

4  1  L  2 


These  represent  useful  information  that  can  be  evaluated  to  validate  the 
candidate  solution. 


b.  Hamiltonian  Evolution  and  Final  Value 

The  Hamiltonian  evolves  in  accordance  with  the  simple  equation, 

H  =  —  (4.12) 

dt 


Previously,  we  dealt  only  with  Hamiltonian  equations  that  had  no  specific 
dependence  on  time  and  therefore  the  time-rate  of  change  was  zero.  In  this 
case  the  magnetic  field  of  the  Earth  introduces  a  time-dependence  into  the 

Hamiltonian.  That  is  Therefore,  the  Hamiltonian  is  not  a  constant  in 

the  interval  under  consideration. 


The  final  value  of  the  Hamiltonian  is  given  by, 


.....  dE  ,de  _ 
H[tf\  +  — +vl  —  =  0 
f  dt,  dt, 


where  the  end  manifold  ( e)  is  written  in  the  standard  form, 
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e{xf,tf)  = 


q,  -sin 


Q2 


V2y 


% 


q4.  -cos 


A  A  \ 


V27 


=  0 


Cfly  +0)oC0S((|) ) 

(£>Zf  —  co0  sin(<|) ) 


Then  by  inspection  we  can  see  that  the  final  value  of  the  Hamiltonian  is  negative 
one. 


H[tf  ]+ 1  =0  — >  H[tf]  =  -1 


Although  the  Hamiltonian  is  not  a  constant  for  the  interval  under  consideration  its 
final  value  represents  a  second  numerical  figure  of  merit  to  validate  the  optimality 
of  the  solution. 


3.  Numerical  Results 

The  numerical  example  for  this  work  was  taken  from  the  Naval 
Postgraduate  School’s  current  small  satellite  program,  “NPSAT  1.”  Designed 
primarily  to  allow  a  hands-on  learning  experience  this  satellite,  still  in  the  design 
phase,  will  provide  three-axis  magnetic  torque  control  with  a  pitch  wheel  for 
increased  stabilization.  The  moment  of  inertias  and  orbital  parameters  used  in 
the  numerical  examples  were  taken  from  NPSAT  1  preliminary  designs.  The 
contribution  of  the  pitch  bias  wheel  will  be  addressed  in  a  later  chapter.  The 
NPSAT  1  data  includes: 


78 


NPSAT 1  Parameters 


Tab 


lx 

5  kg-m2 

ly 

5.1  kg-m2 

lz 

2  kg  m2 

Max  Dipole  Moment 

30  Amp-m2 

Orbital  Altitude 

560  km  (Circular) 

Inclination 

35.4  degrees 

e  5  NPSA  T 1  Parameters  for  Numerical  Simulations 


The  maneuver  selected  for  simulation  is  a  135  degree  roll  (x-axis  slew). 

The  time-optimal  control  solution  for  this  maneuver  demonstrates  a 
surprisingly  clean  bang -bang  structure.  The  solution,  shown  in  Figure  49,  has  10 
control  switches  distributed  among  the  three  axes. 


Figure  49  Time-optimal  Control  Solution  for  Magnetic  Torque  Problem 
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Before  evaluating  the  optimality  of  the  candidate  solution  the 
feasibility  is  evaluated.  The  control  solution  is  propagated  through  a 
separate  ODE  45  dynamics  simulator  to  verify  that  the  candidate  solution 
drives  the  dynamic  system  from  the  initial  state  to  the  final  state.  The 
propagation  results  (Figure  50  &  Figure  51)  show  that  the  control  solution 
does  meet  the  end  point  constraints  and  that  the  estimated  states  closely 
match  those  obtained  during  propagation.  The  original  solution  obtained 
is  shown  in  solid  lines  overlaid  with  the  propagated  states  shown  as  V 
marks  below. 


Ol 


Tim*  (s*el 


Figure  50  Quaternion  Solution  and  Validation  by  Propagation 
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Figure  51  Angular  Rate  Solution  and  Validation  by  Propagation 

The  quaternion  and  angular  rate  histories  are  shown  (Figure  52  &  Figure  53). 
The  maneuver  is  clearly  not  an  eigenaxis  slew.  This  is  evident  from  both  the 
variation  in  the  quaternions  q2  &g3and  the  non-zero  angular  rates  of  co2  &co3. 


Figure  52  Magnetic  Torque  Slew  Time-optimal  Quaternion  History 
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Figure  53  Magnetic  Torque  Slew  Time-optimal  Angular  Rate  History 

Next  we  evaluate  the  optimality  of  the  feasible,  candidate  solution.  The 
switching  functions  are  given  in  equations  (4.1 1).  These  are  plotted  overlaid  with 
the  scaled  control  solution  (Figure  54,  Figure  55,  &  Figure  56). 


Figure  54  Magnetic  Torque  Control  X-axis  Switching  Function  and 

Control  Solution 
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Figure  55  Magnetic  Torque  Control  Y-axis  Switching  Function  and 

Control  Solution 


Figure  56  Magnetic  Torque  Control  Z-axis  Switching  Function  and 

Control  Solution 
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Shown  are  the  switching  functions  (S,)  previously  defined  as  the  partial 

derivative  of  the  Hamiltonian  with  respect  to  the  control  vector.  The  KKT 
multiplier  (Mu)  is  also  shown  plotted  separately  from  the  switching  function.  The 
sum  of  the  switching  function  and  the  KKT  multiplier  is  the  definition  of  the 
minimization  of  the  Hamiltonian,  equation  (4.11),  and  should  be  numerically 
equal  to  zero.  Additionally,  as  before  the  switching  function  and  control  are 
related  be  the  KKT  conditions: 


Uj 


maximum 
\  minimum 
singular 


S(.  <0 

S,>  0 

s,=o 


These  figures  clearly  illustrate  that  the  control  solution  meets  optimality  criteria 
established  by  the  Hamiltonian  minimization. 

The  Hamiltonian  transversality  and  evolution  conditions  are  complicated 
by  the  varying  magnetic  field.  The  Hamiltonian  and  the  predicted  final  value  are 
shown  (Figure  57).  The  numerical  final  value  of  the  Hamiltonian  and  the 
theoretical  final  value  differ  by  only  0.0646.  Therefore  we  conclude  that  the 
candidate  solution  has  met  the  necessary  conditions  for  optimality. 


Figure  57  Magnetic  Torque  Solution  Hamiltonian  and  Final  Value 
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Additionally,  to  validate  our  hypothesis  that  the  time-dependence  of  the 
Hamiltonian  was  due  to  the  varying  magnetic  field  the  problem  was  evaluated 
with  a  constant  magnetic  field  approximation.  The  Hamiltonian  for  the  case  of  a 
constant  magnetic  field  is  shown  in  Figure  58.  The  Hamiltonian  for  this  case  has 
lost  its  apparent  dependence  on  time  and  settled  to  a  value  that  is  numerically 
close  to  the  value  of  -1  that  was  predicted. 


0 

-0.5 

-1 

-1.5 

-2 

o  e 

Constant  B-Field  Hamiltonian 

_ 

_ 

100  200  300  400  500 

Figure  58  Constant  Magnetic  Field  Approximation  Hamiltonian  Solution 
4.  Numerical  Considerations  and  Scaling 

The  linear  scaling  used  throughout  these  algorithms  was  introduced  in  the 
previous  section.  In  this  section  scaling  was  also  added  for  the  Earth’s  magnetic 
field  in  the  form: 

B,  =  kBBi 

The  goal  of  the  scaling  is  to  bring  all  numerical  values  seen  by  the  optimization 
solver  into  the  same  order  of  magnitude.  Scaling  values  were  adjusted  from  an 
unsealed  solution  to  improve  the  quality  of  the  solution  and  then  readjusted  as 
necessary.  Proper  scaling  reduced  computation  time  and  improved  the  accuracy 
of  the  solution.  Additionally,  the  switching  functions  of  properly  scaled  problems 
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behaved  more  closely  to  theoretical  predictions  indicating  that  results  were 
improved  by  proper  scaling. 


Scaling  Effects 

Maneuver  Time  (before 

scaling) 

230.0845  seconds 

Maneuver  Time  (after 

scaling) 

271.1564  seconds 

Error 

15.15% 

Table  6  Effects  of  Scaling  on  Solution  Fidelity 


5.  Conclusions 

In  this  section  the  open  loop  time-optimal  control  for  a  magnetic  torque 
controlled  asymmetric  spacecraft  was  determined.  The  candidate  solution  was 
determined  by  propagation  to  be  a  feasible  solution  to  the  problem.  The 
optimality  of  the  solution  was  validated  through  an  analysis  of  the  Hamiltonian 
minimization,  switching  functions  and  the  behavior  of  the  Hamiltonian.  The 
hitherto  unseen  variation  of  the  Hamiltonian  over  time  was  theorized  to  be 
caused  by  the  time  dependence  of  the  Earth’s  magnetic  field.  When  this 
dependence  was  eliminated  the  Hamiltonian  returned  to  the  constant  values  we 
have  seen  previously.  Therefore,  we  conclude  that  the  solution  is  feasible  and 
meets  the  necessary  conditions  for  optimality. 
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2  Sidi,  M.J.  (1997).  Spacecraft  Dynamics  and  Control.  Cambridge  University  Press. 
New  York,  NY.. 

3  Presented  in  the  class  notes  of  Professor  B.  Leonard  AA-3818  /  Fall  ’04.  Spacecraft 
Dynamics  and  Control.  US  Naval  Postgraduate  School,  Monterey,  CA. 
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V.  NPSAT  1  CONTROL  SYSTEM 


A.  INTRODUCTION 

NPSAT  1  is  a  small  satellite  design  project  currently  in  work  at  the  U.S. 
Naval  Postgraduate  School  (Figure  59).  It  was  conceived  as  a  three-axis 
stabilized  magnetic  torque  controlled  satellite.  Later  in  the  design  process,  a 
pitch  bias  wheel  was  added  to  improve  stability  and  reliability.  In  this  section  we 
explore  the  time-optimal  reorientation  of  this  small  asymmetric  satellite  where 
both  the  three  torque  rods  and  the  pitch  wheel  are  available  as  torque  generating 
devices. 


Figure  59  NPSA  T 1  Conceptual  Image  courtesy  of  Dan  Sakoda 
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B.  PROBLEM  FORMULATION 


In  this  case  the  spacecraft  state  is  defined  by  its  position,  angular  velocity 
and  the  angular  momentum  of  the  pitch  wheel. 


x  = 


Q 

co 

h 


The  position  is  represented  by  a  four-element  quaternion  vector;  we  have 
previously  adopted  the  convention  that  the  fourth  element  of  the  quaternion 
vector  is  the  scalar  quantity.  The  angular  velocity  is  in  body  coordinates  with 
respect  to  the  Newtonian  frame.  The  pitch  wheel  is  assumed  to  be  aligned  with 
the  spacecraft  #2  principal  axis.  This  assumption  simplifies  the  formulation  of  the 
rotational  dynamics  equations.  The  orbital  parameters  are  as  given  in  Table  5. 


The  quaternion  kinematics  equations  are  well  known  and  shown  in 
equation  (3.33).  These  equations  are  unchanged  by  the  addition  of  the  pitch  bias 
wheel.  The  rotational  dynamics  however,  are  now  the  result  of  the  external 
torque  generated  by  the  interaction  of  the  magnetic  field  and  the  torque  rods  in 
addition  to  the  pitch  wheel  angular  momentum  and  torque. 

The  time  rate  of  change  of  angular  momentum  can  be  expressed  as1, 


d 

d  n 

—H , 

— 

—Hq 

L dt  s\ 

N 

[dt  s\ 

+  (5  X  Hr 


(5.1) 


where  Hs  is  the  total  angular  momentum  of  the  spacecraft-wheel  system  and  is 

expressed  in  the  body  frame.  Then  by  assuming  that  the  wheel’s  center  of  mass 
is  collocated  with  the  spacecraft  center  of  mass  we  can  express  the  system 
angular  momentum  as, 

H  =  lB%b  +  lw%w  (5.2) 


where  lB&lw  are  the  body  and  wheel  moments  of  inertia  respectively,  and  the 
angular  velocity  of  the  wheel  with  respect  to  the  Newtonian  frame  is  given  by, 


Then  the  angular  momentum  of  the  system  can  be  written  as, 


By  defining, 


Hq 


lBN(Db 


+  LN<»b 


+lwbtiw 


J  —  I B  +  l\W 


we  have, 

Hs  =  JN(5b  +  lwb65w 


(5.3) 


This  result  is  easily  derived  under  the  assumptions  stated.  Kane2  provides  a 
detailed  derivation  to  show  that  this  result  holds  for  any  configuration. 

Then,  referring  to  equation  (5.1)  we  have, 

Mext  =  JNcbb  +  lw  bcbw  +  N(db  x(JN<5b+lw  b(5w) 


By  defining, 


•w 


b(DW 


we  can  write , 

Mext  =  JN (5b  +  hw  +  N(£>bxJN(5b  +  hw 


If  we  allow  Mext  to  be  the  sum  of  disturbance  torques  and  the  torque  generated 
by  the  interaction  of  the  Earth’s  magnetic  field  and  the  magnetic  torque  rods, 
setting  the  disturbance  torque  to  zero  we  can  write, 


Jco  +(3x  Jco  =  mx  B-hw  ~(5xhw 
where, 


1 

'  O 

_ i 

1 

o 

1 _ 

hw  — 

hw 

>  hw  - 

hw 

l 

O 

i _ 

1 

o 

_ 1 

and  co  =  wcob 


Carrying  out  the  cross  products  and  rotational  transformations  previously  defined 
gives  the  following  result: 
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JiCOx  -  ( J2  -  J3  )coycoz  =  m2B3  -  m3B2  +  co  zhw 
J2coy  -  (J3  -J,  )co/oz  =  m3B ;  -  m,B3  -  hw 
J3coz  -  ( J,  -  J2  )coxcoy  =  m,B2  -  m2B:  -  co xhw 

where  the  Earth’s  Magnetic  Field  (6,,%^)  is  in  the  spacecraft  body  frame  and 
was  previously  defined  by  equation  (4.6).  The  four-element  control  vector  is 
chosen  as: 

~m,~ 

m7 

u=  2  (5.4) 

m3 

_uw_ 

Box  constraints  are  imposed  on  the  control  elements  based  on  the  physical 
limitations  of  the  control  elements  selected.  The  components  under 
consideration  in  this  model  have  the  physical  characteristics  shown  (Table  7). 


Component  Characteristics 

Torque  Rods 

Maximum  Magnetic 

Dipole  Moment 

30  Amp*m2 

Pitch  Wheel 

Maximum  Angular 

Momentum 

18  N*m*sec 

Pitch  Wheel 

Maximum  Torque 

30  mN*m 

Table  7  NPSAT 1  Simulation  Component  Characteristics 


Pitch  wheel  angular  momentum  limits  are  imposed  as  a  state  constraint.  Other 
values  are  imposed  as  control  constraints. 

C.  TIME-OPTIMAL  NPSAT  1  SLEW 

In  this  section  we  consider  the  time-optimal  reorientation  of  an  asymmetric 
spacecraft  controlled  by  a  combination  of  magnetic  torque  rods  and  a  pitch 
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momentum  wheel.  The  maneuver  is  defined  as  rest-to-rest  in  the  orbit  frame 
where  the  initial  attitude  and  attitude  rates  are  known.  The  final  attitude  and 
attitude  rates  are  determined  by  the  eigenaxis  of  rotation  and  the  rotation  angle. 
The  angular  momentum  of  the  wheel  is  the  final  state  variable.  For  this  variable 
and  control  combination  there  are  four  possible  options.  First,  the  initial  wheel 
angular  momentum  may  be  left  unspecified  as  an  optimization  variable  to  be 
determined.  Then,  if  wheel  torque  does  not  equal  zero,  the  algorithm  selects  the 
initial  wheel  speed  and  control  history  for  the  minimum  time  maneuver.  This  is 
the  case  that  is  numerically  evaluated  below.  The  wheel  speed  can  also  be  fixed 
at  the  end  points.  In  this  case  a  non- zero  toque  limit  will  achieve  the  time  optimal 
solution  within  the  constraints  provided.  Under  the  condition  of  zero  torque, 
wheel  speed  free,  the  algorithm  will  determine  the  optimal,  constant  wheel  speed 
for  the  desired  maneuver.  Finally,  if  wheel  speed  and  torque  are  set  to  zero  the 
results  match  those  previously  obtained  for  the  magnetic  torque  control  section. 
It  is  important  to  note  that  these  changes  in  boundary  conditions  do  not  affect  our 
ability  to  minimize  the  objective  function  and  obtain  valid  solutions. 

1.  Problem  Statement 

Determine  the  controls  [u* ,u2* ,uz* ,uA*]  that  drive  the  spacecraft  from  its 
initial  rest  position,  given  by  [xq] to  its  final  rest  position  given  by  [xf] while 
minimizing  the  cost  function: 


where  we  have  used  the  Mayer  form  of  the  cost  function,  subject  to  the  following 
constraints: 

Control  /  State  Constraint:  The  control  and  state  constraints  are  defined  in 
the  standard  form, 

hL{t)<h{u,t)<hu{t) 
then  our  control  constraints  are  written, 
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-30  <  Uj  <30  (Amp-m2)  /  =  1,2,3 
-30x1 0  3  <  u4  <30x1 0“3  (N  m) 

and  the  state  constraint  is  given  by: 

-0.18  <x8  <0.18  (N-m  sec) 


Dynamic  Constraints: 


=^N4-^3+»3(?2] 

4  =  -  co3Qi  ] 

<73=^[-°V72  +“2^1  +©3<74] 

«x  =  -7- (m2e3  - m3B2  -K>zhw)  +  kpyaz, 

Jx 

«V  =  -)-( m3B,  ~  mv B3  +  K )  +  M#,. 

Jy 

=  J-(^e2  -  m2B ;  +  toA)+/wv 


Jy~JZ 

J, 


k  _JZ~JX 
2~ 

, 


(5.5) 


hw  =  u* 
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2.  Solving  the  Optimal  Control  Problem 

The  Hamiltonian  is  given  by, 


X  X 

H  =  y  (®xq4  -ayq3  +coz9,2  +“oQ3  )  +  y  (coxg3  -noy  -co zq,  +co 0q, )  + 
X  X 

y(-®*92  +“/7i  +“2Qz-w0Qi)  +  yHv?i  -<V& -“zQs  -®0<72)  + 


1 


1 


-T  (U2e3  -  U3B2  -®z^  )  +  *W°Z  +  V  T  (U3ei  -  UA  + W4)  +  M^z  + 


J 


V \-(uA-u2B,+03xhw)  +  k3 C0x0)y  [  +  \ u4 


J. 


The  subscripts  on  the  Lagrange  multipliers  have  been  chosen  for  bookkeeping 
purposes.  The  control  vector  is  defined  in  equation  (5.4). 


Following  the  solution  method  previously  established  we  begin  by 
minimizing  the  Hamiltonian  subject  to  the  control  constraint  set.  Then  the 
necessary  conditions  for  the  Hamiltonian  minimization  are  as  follows: 


X 


-^b-^b3  +  M  =0 


J 


J 


Jx  Jz 

-j~B,-^B2  +  p.3=0 

y  u  x 


K 

j , 


~+Xh  +  m  -  o 


muj 


(5.6) 


These  relationships  will  be  evaluated  to  validate  the  candidate  solution. 

We  saw  that  in  the  case  of  the  magnetic  torque  problem  the  magnetic  field 
of  the  Earth  introduced  a  time-dependence  into  the  Hamiltonian.  Therefore,  the 
Hamiltonian  was  not  a  constant  in  the  interval  under  consideration.  For  the 
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current  problem  the  Earth’s  magnetic  field  is  still  a  factor  in  the  Hamiltonian. 
However,  in  this  case,  we  will  see  that  the  effect  of  the  magnetic  field  is  reduced 
by  the  introduction  of  the  pitch  wheel,  and  the  time  rate  of  change  of  the 
Hamiltonian  is  reduced  significantly. 

The  final  value  of  the  Hamiltonian  is  again  given  by, 


.......  dE  tde  _ 

H[tf]  +  —  +  vl  —  =0 


at 


dtf 


where  the  end  manifold  ( e)  is  written  in  the  standard  form  previously  established, 

e{xf,tf  )=0 

Then  by  inspection  we  can  see  that  the  final  value  of  the  Hamiltonian  is  negative 
one. 

H[f,]+1=0-»H[f,]  =  -1 


3.  Numerical  Results 

The  numerical  example  for  this  work  was  again  taken  from  the  Naval 
Postgraduate  School’s  current  small  satellite  program,  “NPSAT  1.”  The  moment 
of  inertias  and  orbital  parameters  were  previously  established  in  Table  5  and  are 
assumed  to  be  the  system  moments  of  inertia.  The  control  parameters  are  given 
in  Table  7. 

The  maneuver  selected  for  simulation  is  a  135  degree  roll  (x-axis  slew) 
and  is  a  rest-to-rest  maneuver  in  the  orbit  frame.  The  initial  and  final  pitch  wheel 
angular  momentums  are  free  within  the  state  bounds  as  an  optimization 
parameter.  The  time-optimal  control  solution  is  shown  in  Figure  60  and  Figure 
61. 
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Mag-Wheel  Magnetic  Dipole  Moment  Control  Solution 
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Figure  60  NPSAT 1  Torque  Rod  &  Wheel  Time-optimal  Control  Solution 

(V 


Figure  61  NPSAT  1  Pitch  Wheel  Torque  Time-optimal  Control  Solution  (2) 


There  are  several  interesting  characteristics  to  this  solution  not  the  least  of 
which  is  an  apparent  singular  arc  in  the  switching  structure  for  the  pitch  wheel 
torque  ( u4).  Before  evaluating  the  optimality  of  the  candidate  solution  the 

feasibility  is  evaluated.  The  control  solution  is  again  propagated  through  a 
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separate  ODE  45  dynamics  simulator  to  verify  that  the  candidate  solution  drives 
the  dynamic  system  from  the  initial  state  to  the  final  state.  Previous  control 
solutions  were  propagated  using  a  linear  interpolation.  In  this  case  a  piecewise 
cubic  hermite  interpolating  polynomial  (pchip)  produced  more  accurate  results. 
The  propagation  results  figure  62,  Figure  63,  and  Figure  64)  show  that  the 
control  solution  does  meet  the  end  point  constraints  and  that  the  estimated  states 
closely  match  those  obtained  during  propagation.  The  original  solution  obtained 
is  shown  in  solid  lines  overlaid  with  the  propagated  states  shown  as  V  marks 
below.  Therefore  we  conclude  that  the  control  solution  is  feasible. 

The  state  histories  are  shown  figure  65,  Figure  66,  and  Figure  67).  It 
comes  as  no  surprise  that  the  time-optimal  maneuver  is  not  an  eigenaxis 
maneuver.  This  is  evident  from  both  the  variation  in  the  quaternions  q2  &g3and 
the  non-zero  angular  rates  of  co2  &co3. 
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Figure  63  Angular  Rate  Solution  and  Validation  by  Propagation 


Figure  64  Pitch  Wheel  Momentum  Solution  and  Validation  by 

Propagation 
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Figure  66  NPSA  T 1  Slew  Time-optimal  Angular  Rate  History 
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Figure  67  NPSAT 1  Slew  Time-optimal  Pitch  Wheel  Angular  Momentum 

History 

Next  we  evaluate  the  optimality  of  the  feasible,  candidate  solution.  The 
switching  functions  are  given  in  equations  (4.1 1).  These  are  plotted  overlaid  with 
the  scaled  control  solution  (Figure  68,  Figure  69,  Figure  70,  and  Figure  71 ). 


Figure  68  NPSAT  1  Control  Switching  Function  and  Control  Solution  (1) 
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Figure  69  NPSAT 1  Control  Switching  Function  and  Control  Solution  (2) 
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Figure  70  NPSAT  1  Control  Switching  Function  and  Control  Solution  (3) 
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Figure  71  NPSAT 1  Control  Switching  Function  and  Control  Solution  (4) 


Shown  are  the  switching  functions  (S,)  previously  defined  as  the  partial 

derivative  of  the  Hamiltonian  with  respect  to  the  control  vector.  The  KKT 
multiplier  (Mu)  is  also  shown  plotted  separately  from  the  switching  function.  The 
sum  of  the  switching  function  and  the  KKT  multiplier  is  the  definition  of  the 
minimization  of  the  Hamiltonian,  equation  (4.11),  and  should  be  numerically 
equal  to  zero.  As  before  the  switching  function  and  control  are  related  by  the 
KKT  conditions: 


u. 


maximum 
\  minimum 
singular 


S(.  <0 

Si  >0 

Sj  -0 


The  Hamiltonian  for  this  problem  is  shown  in  Figure  72.  As  we  alluded  to 
earlier  the  time  dependence  observed  in  the  magnetic  torque  control  problem 
appears  to  have  been  eliminated.  In  fact,  it  has  been  reduced  to  the  extent  that  it 
is  no  longer  visible.  By  constraining  the  pitch  wheel  to  zero  torque  and  zero 
angular  momentum  the  previous  magnetic  torque  results  are  obtained. 
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Hamiltonian  Evolution 


Figure  72  NPSAT 1  Torque  Rod  &  Wheel  Problem  Hamiltonian  Evolution 

These  figures  clearly  illustrate  that  the  control  solution  meets  optimality 
criteria  established  by  the  Hamiltonian  minimization.  A  closer  look  at  the 
Hamiltonian  minimization  with  respect  to  the  pitch  wheel  torque  is  given  in  Figure 
73.  The  switching  function  appears  singular  during  the  period  in  which  the  pitch 
wheel  torque  is  zero. 


Figure  73  NPSAT  1  Control  Switching  Function  for  Wheel  Torque  (Close 

up) 
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4.  Numerical  Considerations  and  Scaling 

We  continue  with  the  linear  scaling  established  previously.  In  this  section 
scaling  was  also  added  for  the  pitch  wheel  parameters  in  the  form: 

h=khh 

h=^h 

kt 

As  before  the  goal  of  the  scaling  is  to  bring  all  numerical  values  seen  by  the 
optimization  solver  into  the  same  order  of  magnitude.  Scaling  values  were 
adjusted  from  an  unsealed  solution  to  improve  the  quality  of  the  solution  and  then 
readjusted  as  necessary. 

5.  Conclusions 

In  this  section  the  open  loop  time-optimal  control  for  an  asymmetric 
spacecraft  equipped  with  three  magnetic  torque  rods  and  a  pitch  wheel  was 
determined.  The  candidate  solution  was  determined  by  propagation  to  be  a 
feasible  solution  to  the  problem.  The  optimality  of  the  solution  was  validated 
through  an  analysis  of  the  Hamiltonian  minimization,  switching  functions  and  the 
behavior  of  the  Hamiltonian.  Therefore,  we  conclude  that  the  solution  is  feasible 
and  meets  the  necessary  conditions  for  optimality. 

By  allowing  the  use  of  the  pitch  wheel  as  a  spacecraft  control,  the  time 
required  for  the  reorientation  maneuver  was  significantly  reduced.  In  the  original 
configuration,  with  only  magnetic  torque  control  available,  the  time  required  for 
the  optimal  maneuver  was  271 .1564  seconds.  In  the  new  configuration  the  same 
maneuver  required  only  141.9933  seconds.  This  represents  a  47.6  percent 
reduction  in  time  required  for  the  same  maneuver. 
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ENDNOTES 


1  Wie,  B.  (1998).  Space  Vehicle  Dynamics  and  Control.  AIAA  Education  Series, 
American  Institute  of  Aeronautics  and  Astronautics,  Inc.  Reston,  VA. 

2  Kane,  T.  (1983).  Spacecraft  Dynamics.  McGraw-Hill  Inc.  New  York,  NY. 
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VI.  CONTROL  MOMENT  GYRO  CONTROL  SYSTEMS 


A.  INTRODUCTION 

Control  moment  gyros  (CMGs)  are  well  known  for  their  large  torque 
generating  capability.  However,  in  the  past,  because  of  their  mechanical 
complexity,  cost  and  weight  they  have  been  restricted  to  large  satellites  with  high 
agility  requirements.  Recent  research  in  the  design  of  smaller,  less  expensive 
CMGs  has  created  a  renewed  interest  in  CMG  spacecraft  control.  However,  the 
primary  focus  of  this  large  body  of  work  has  been  singularity  avoidance12  3  4.  In 
this  section  we  calculate  the  time-optimal  solution  for  an  asymmetric  spacecraft 
using  the  much  studied  four-CMG  pyramid  configuration.  In  addition  to 
completing  a  time-optimal  maneuver  we  will  show  that  the  problem  formulation 
generates  a  singularity  free  solution. 

B.  CONTROL  MOMENT  GYRO  BASICS 

A  control  moment  gyro  contains  a  flywheel  which  spins  at  a  constant  rate. 
The  spin  axis  of  the  flywheel  is  connected  to  a  gimbal  which  allows  reorientation 
of  the  spin  axis.  Therefore,  the  direction  of  the  angular  momentum  vector  can  be 
changed  with  respect  to  the  spacecraft  body  frame  by  gimballing  (Figure  74). 


Single-gimbal  Control  Moment  Gryo 


Torque  Vector 


Gimbal  M 


Gyro  Motor 


Flywheel 


Angular  Momentum  Vector 


Figure  74  Single-Gimbal  Control  Moment  Gyro  (After  Ref. [5]) 
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For  single-gimbal  control  moment  gyros  (SGCMG),  the  spin  axis  reorientation  is 
restricted  to  a  plane  perpendicular  to  the  gimbal  axis.  The  advantage  of  CMGs  is 
that  small  gimballing  results  in  a  large  control  torque  on  the  spacecraft.  This 
resultant  torque  is  orthogonal  to  both  the  gimbal  and  spin  axes.6  The  torque  from 
the  CMG  is  given  by: 

^cmg  =  hx§  (6.1 ) 

where  h  is  the  angular  momentum  vector  and  has  units  of  Newton-meter- 

seconds,  5*  is  the  gimbal  angle  rate  with  units  of  radians  per  second.  This 
torque  amplification  property  makes  CMGs  desirable  for  applications  that  require 
spacecraft  agility. 


C.  PROBLEM  FORMULATION 

Following  a  similar  development  to  that  used  in  Chapter  V,  recall  that  the 
time  rate  of  change  of  angular  momentum  can  be  expressed  as, 


Mext  = 


d 

d 

—  He 

— 

[  dt  s\ 

N 

[dt  s\ 

+  (5  X  Hr 


(6.2) 


where  Hs  is  the  total  angular  momentum  of  the  spacecrat-CMG  system  and  is 

expressed  in  the  spacecraft  body  frame.  Assuming  that  the  CMG  center  of  mass 
is  collocated  with  the  spacecraft  center  of  mass  we  can  express  the  system 
angular  momentum  as: 

h  -  iB  (£>  +  iCMG  co 


where  lB  and  lCMG  are  the  spacecraft  and  CMG  system  moment  of  inertias 

respectively  and  the  angular  velocity  of  the  CMG  with  respect  to  the  Newtonian 
frame  is  given  by, 

n65cmg  =  wcob  +  b(3CMG 

Then  the  angular  momentum  of  the  system  can  be  written  as, 
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Once  again,  by  defining, 


J  —  Ib  +  I CMG 


we  have, 

Hs  =  J%b  +  lCMGbc »CMG 

Then  referring  to  equation  (6.2)  we  have, 

Mext  =  JNd5b  +  lCMGb<bCMG+  %bx(JN&b  +  lCMGb<aCMG) 

Finally,  by  defining, 

—  I  b(Ti  I-, 

~  CMG  ^  _  '  CMG 

we  obtain  the  following: 

Mext  =J%b  +  hCMG+  %bx  J%b  +  hCMG 

If  we  allow  wcofa  =  co  and  Mext  to  be  all  disturbance  torque,  then  setting 
disturbance  torque  equal  to  zero  we  write, 

J(5  +  co  x  J(5  =  -hCMG  -  63  x  hCMG  (6.3) 

The  CMG  angular  momentum  vector  (/?)*  is  a  function  of  the  gimbal 
angles  (5  )  in  the  body  frame  and  for  multiple  CMG  configurations,  a  function  of 
the  configuration.  Here  we  consider  the  four-CMG  pyramid  configuration  (Figure 
75).  Each  face  of  the  pyramid  is  inclined  from  the  horizontal  by  a  skew  angle 
(P  ).  The  four  CMGs  have  gimbal  axes  orthogonal  to  the  pyramid  faces  and  so 

are  constrained  to  gimbal  on  the  faces  of  the  pyramid.  We  have  selected  a  skew 
angle,  p  =  54.73  degrees.  For  CMGs  with  equal  angular  momentum  about  the 


d_r 

^  nCMG 


*  We  have  adopted  the  notation  lower  case  ‘h’  as  CMG  angular  momentum  and  subsequently 
dropped  the  subscript. 


spin  axis  this  configuration  results  in  a  nearly  spherical  momentum  envelope. 
Additionally,  this  configuration  has  been  extensively  studied  in  the  literature.7 


Reference 

Frame 


CMG  #2 
Gimbal  Axis 

S, 

Pitch  Axis 

CMG  #1 
Gimbal 
Axis 


Figure  75  Pyramid  Mounting  Arrangement  of  Four  Single  Gimbal  CMGs 

(from  Ref. [8]) 

Then  referencing  Figure  75,  we  can  write  the  total  CMG  angular 
momentum  expressed  in  the  body  frame  as: 

-cpsinSJ  [  -cosS2  1  [cpsinSgl  cosS4 

h=  cos81  +  -c(3sin82  +  -cos83  +  cpsinS4  (6.4) 

spsinS^  sp  sinS2  J  |_sp  sin S3  J  |_sp  sinS4 

where  cf>  =  cosp  and  sf>  =  sin  p  .  The  angular  momentum  magnitude  is  set  to 
unity  without  loss  of  generality  {h0=  1).  Then  the  time  of  rate  of  change  of 
angular  momentum  may  be  written  as: 

-cpcosS1  sinS2  cpcosS3 

— sinS!  -cp  cos82  sin83 

sp  cosS!  spcos82  spcos83 

A 
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The  matrix  A  is  defined  from  equation  (6.4)  as, 


and  is  a  3 xn  Jacobian  matrix  where  n  is  the  number  of  CMGs  in  the 
configuration.  Expanding  equation  (6.3)  showst, 

~  K  -  ^z)“y«z  =  -®2  h,  +co3 1\ 

Jy(bV  -  ( Jz  -J>Pz  =  ~K  +  ®A 

JPZ  ~(JX  ~Jy  K0)y  =  -/?3  -  0)1  ^  +  0)2  ^ 

where  the  values  of  fyare  established  in  equation  (6.4)  and  the  values  of  h:  are 

established  in  equation  (6.5).  A  lengthy  direct  substitution  allows  us  to  write  the 
time  rate  of  change  of  angular  velocity  in  the  standard  form: 


where  the  state  vector  and  control  vector  are  defined  as: 

state:  [q,co,6]r  e  M11 
control:  [u]r  e  K4 

D.  CMG  TIME-OPTIMAL  SLEW  MANEUVERS 

In  this  section  we  consider  the  time-optimal  reorientation  of  an  asymmetric 
spacecraft  controlled  by  control  moment  gyros.  In  order  to  simplify  the 
formulation  the  spacecraft  is  assumed  to  be  in  inertial  space.  The  maneuver  is 
defined  as  rest-to-rest  in  the  inertial  frame  where  the  initial  attitude  and  attitude 
rates  are  known.  The  final  attitude  and  attitude  rates  are  specified.  The  gimbal 
angle  of  the  CMG,  the  final  state  variable  is  left  free  as  an  optimization  variable. 


t  This  derivation  is  also  found  in  reference  [6]  where  it  is  used  in  the  development  of  a 
feedback  control  law. 
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1.  Problem  Statement 

Determine  the  controls  [u^,u2*,u3*,u4*]  that  drive  the  spacecraft  from  its 
initial  rest  position,  given  by  [xq] to  its  final  rest  position  given  by  [xf] while 
minimizing  the  cost  function: 

J(x(-),u(-),tf)  =  tf-t0 

where  we  have  used  the  Mayer  form  of  the  cost  function,  subject  to  the  following 
constraints: 

Control  /  State  Constraint:  The  control  and  state  constraints  are  defined  in 
the  standard  form, 

hL{t)<h{u,t)<hu{t) 
then  our  control  constraints  are  written, 

-1  <Uj<  1  rad/sec  /  =  1,2, 3, 4 

and  the  state  constraint  is  given  by: 

-7i  <  S(.  <  71  rad  /  =  1 ,2,3,4 

Physically,  these  equate  to  a  CMG  that  is  capable  of  full  360  degree  rotation 
about  its  gimbal  axis  at  a  rate  of  1  radian  per  second.  With  a  unit  angular 
momentum  then  from  equation  (6.1 )  the  maximum  torque  output  from  the  CMG  is 
1  Newton-meter. 
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Dynamic  Constraints: 


q  =  —£lq 
-  2  =- 

Jy  _  Jz  h-i 

cox=^— coyco2--^ 

J z  —  Jx  ho 

®r  =  ^p®A-j 


_ho 

J ,  ~'x”y  J, 


CO,  =— - — COG)  —  — 


tK^-coA) 

Jx 

— (C0z/7i-C0x/73) 
Jy 

— (cox/?2-coy/?i) 


(6.6) 


8  /  =  U;  i=  1 . 4 

where  the  quaternion  kinematics  equations  shown  have  been  previously  defined. 
The  quantity  /?,  is  defined  in  equation  (6.5)  in  terms  of  the  state  6  and  the  control 

and  the  quantity  /?,  is  defined  in  equation  (6.4)  in  terms  of  the  state  5  . 


2.  Solving  the  Optimal  Control  Problem 

As  before,  the  first  step  in  solving  the  optimal  control  problem  is  to  form 
the  Hamiltonian.  Basic  format  of  the  Hamiltonian  is  repeated  here*, 

H(k,  x,  u,  t)  =  F(x,  u,  t)  +  XTf(x,  u,  t) 

Since  the  cost  functional  was  formulated  without  a  Lagrange  cost  term  the 
Hamiltonian  reduces  to  the  following. 

H{\,x,u,t)  =  XTf{x,u,t)  (6.7) 

Substituting  equation  (6.6)  into  equation  (6.7)  gives  the  Hamiltonian  for  the  CMG 
spacecraft  system. 


i  Recall  that  in  this  notation  F  is  the  Lagrange  (running)  cost  and  f  are  the  state  dynamics. 
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"  =  W 


-Clq 
o  =— 

v.  y 

y 


+  X 


Jy-J> 


V 


“yroz-y-y  K^-®A) 


+ 


^  I  ^7 -(OA-y-^r®^) 

Jy  Jy  Jy 

f- Jx  .  Jy  (D/Oy -COy/?i) 


V 


V  J. 


4  4 


+ 


+ 


(6.8) 


y 


Vu 


Recall  from  equation  (6.7)  that  the  Hamiltonian  is  written  in  terms  of  the 
Lagrange  multipliers,  state  vector,  control  vector  and  time.  In  order  for  equation 
(6.8)  to  be  rigorously  correct  substitutions  from  equation  (6.4)  and  equation  (6.5) 
are  required.  However,  this  lengthy  substitution  is  omitted. 

Since  the  control  is  constrained,  Hamiltonian  minimization  is  accomplished 
by  adjoining  the  control  constraint  equations  to  the  Hamiltonian  as, 

H{X,,x,u,t,\x)  =  H(k, x,  q  t)  +  \iTh{u,  t)  (6.9) 


Then,  by  differentiation  of  equation  (6.8),  with  respect  to  the  control  vector  and 
substituting  from  equation  (6.4)  and  (6.5),  we  can  write  the  necessary  conditions 
for  the  Hamiltonian  minimization  as  follows. 


—  =— cos(3  cos81  +  -^sin81  sin  (3  cosS^A,^  +p8i  =0 
dU  1  Jx  Jy  Jz 


dH 

du2 

dH 

du^ 


k  sin52  +  ^cos  p  C0S52  -  ksin  p  CosS2  +  ^  +  -  0 


J, 


-cosp  cosS3 


S 

j 


■sinSq 


K 

j. 


■sin(3  cos83  +  \  +(i  =0 


-^-  =  —  sin84 --^-cospcos84 --^sin  pcos84  +  Ag4  +pS4  =0 
duA  Jv  Jw  i 


J ' 


MU: 


(6.10) 


These  conditions  will  be  evaluated  to  validate  the  candidate  solution. 
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Inspection  of  the  Hamiltonian  equation  reveals  no  direct  dependence  on 
time.  Therefore  the  time  rate  of  change  of  the  Hamiltonian  is  zero  and  we  expect 
a  constant  valued  Hamiltonian  for  the  interval  under  consideration. 

The  final  value  of  the  Hamiltonian  is  given  by, 

dE  tde 

H[tf  ]  +  — +1/  —  =  0 
'  dtf  dtf 

where  the  end  manifold  ( e)  is  written  in  the  standard  form  previously  established, 

e{xf,tf  )=0 

Then  by  inspection  we  can  see  that  the  final  value  of  the  Hamiltonian  is  negative 
one. 


H[tf]+J\=0^H[tf]  =  -J\  (6.11) 

Thus  we  expect  a  Hamiltonian  with  a  constant  value  of  -1  over  the  interval  of  the 
maneuver. 


3.  Numerical  Results 

For  this  numerical  example  we  used  the  asymmetric  moment  of  inertia  of 
the  NPSAT  1  small  satellite  design  previously  established  (Table  5)  and 
assumed  to  be  the  system  moments  of  inertia.  The  orbital  parameters  were 
neglected  since  the  model  was  assumed  in  inertial  space.  The  physical 
parameters  for  the  CMG  are  summarized  in  Table  8. 


NPSAT  1  SIMULATED  CMG  SYSTEM 

Gimbal  Rotation  Range 

360  deg 

Maximum  Gimbal  Rotation  Rate 

1  rad/sec 

Maximum  Angular  Momentum 

1  N-m-sec 

Maximum  Output  Torque 

1  N-m 

Table  8  Simulated  CMG  Characteristics 
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The  maneuver  selected  for  simulation  is  a  135  degree  roll  (x-axis  slew) 
and  is  a  rest-to-rest  maneuver  in  the  inertial  frame.  The  initial  and  final  CMG 
gimbal  angles  are  free  within  the  state  bounds.  Therefore,  we  can  write  the  end 
point  conditions  as  follows: 


x. 


x  =  i  q >  q2 ,  q3 q4  -“i  >«*  ^  §3 §4  ]7 r 

x0  =  [0,0, 0,1 ,0,0, 0,8^  ,8%  \  ,84o  f 


sin 


135 


,0,0, cos 


135 


-iT 


,0,0,0, 81  ,S2  ,S3  ,S4 


The  time-optimal  control  solution  is  shown  in  Figure  76.  Our  engineering 
intuition  would  lead  us  to  expect  a  bang-bang  solution  for  a  time-optimal 
maneuver.  The  control  solution  indicates  that  the  solution  is  a  bang-bang 
response  from  all  four  CMGs  with  5  control  switches.  This  is  a  surprisingly  clean 
solution  considering  the  complexity  of  the  CMG  dynamics  derived  above. 

As  we  have  seen  previously,  the  feasibility  of  the  candidate  solution  is 
evaluated  propagating  the  control  solution  through  a  separate  ODE  45  dynamics 
simulator  to  verify  that  the  candidate  solution  drives  the  dynamic  system  from  the 
initial  state  to  the  final  state.  A  linear  interpolation  was  used  in  the  propagation 
sub-routine.  The  propagation  results  (Figure  77,  Figure  78,  and  Figure  79)  show 
that  the  control  solution  does  meet  the  end  point  constraints  and  that  the 
estimated  states  closely  match  those  obtained  during  propagation.  The  original 
solution  obtained  is  shown  in  solid  lines  overlaid  with  the  propagated  states 
shown  as  V  marks  below.  Therefore  we  conclude  that  the  control  solution  is 
feasible. 
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CMG  Control  History 


Time  (sec) 


Figure  76  CMG  Time-optimal  Control  Solution 


Figure  77  CMG  Solution  Quaternion  History  Validation  by  Propagation 
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Figure  78  CMG  Solution  Angular  Rate  History  Validation  by  Propagation 


Figure  79  CMG  Solution  Gimbal  Angle  History  Validation  by  Propagation 
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The  state  histories  are  shown  (Figure  80,  Figure  81 ,  and  Figure  82).  Once 
again  we  see  that  the  time-optimal  maneuver  is  not  an  eigenaxis  maneuver.  This 
is  evident  from  both  the  variation  in  the  quaternions  g2&g3and  the  non-zero 

angular  rates  of  co2  &co3.  Additionally,  the  gimbal  angle  history  demonstrates 


Figure  81  CMG  Slew  Time-optimal  Angular  Rate  History 
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the  state  is  reasonably  well  behaved  and  within  the  state  constraints  [-71,71] 

imposed  in  the  problem  formulation.  The  time  required  to  complete  the 
maneuver  is  4.24  seconds.  This,  as  expected,  is  faster  than  the  previous 
idealized  actuator  solution  for  this  moment  of  inertia  configuration  (see  Chapter 
III).  The  torque  available  from  the  four  actuators  is  higher  than  the  previous 
idealized  actuator  toque  numerical  example  and  therefore  we  expect  faster 
maneuvering.  So  our  engineering  judgment  leads  us  to  believe  that  the  solution 
is  correct. 

Next  we  evaluate  the  optimality  of  the  feasible,  candidate  solution.  The 
switching  functions  are  given  in  equations  (6.10).  These  are  plotted  overlaid  with 
the  control  solution  (Figure  83  through  Figure  86).  Shown  are  the  switching 
functions  (S,  )  previously  defined  as  the  partial  derivative  of  the  Hamiltonian  with 

respect  to  the  control  vector.  The  KKT  multiplier  (Mu)  is  also  shown  plotted 
separately  from  the  switching  function.  The  sum  of  the  switching  function  and 
the  KKT  multiplier  is  the  definition  of  the  minimization  of  the  Hamiltonian, 
equation(4.10),  and  should  be  numerically  equal  to  zero. 
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Figure  83  CMG  Control  Switching  Function  and  Control  Solution 


Figure  84  CMG  Control  Switching  Function  and  Control  Solution  (2) 
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Figure  85  CMG  Control  Switching  Function  and  Control  Solution  (3) 


Figure  86  CMG  Control  Switching  Function  and  Control  Solution  (4) 
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As  before  the  switching  function  and  control  are  related  by  the  KKT  conditions: 


u. 


maximum 
•  minimum 
singular 


Si  <0 

Si  >0 

Sj  -0 


These  figures  clearly  illustrate  that  the  control  solution  meets  the  necessary 
conditions  for  optimality  established  by  the  Hamiltonian  minimization. 

The  Hamiltonian  transversality  and  evolution  conditions  were  defined  in 
equations  (4.12)  and  (6.1 1).  The  computed  Hamiltonian  is  shown  in  Figure  87.  It 
is  clear  that  the  Hamiltonian  is  numerically  well  behaved  and  meets  the 
necessary  conditions  for  optimality. 


Figure  87  CMG  Time-optimal  Slew  Solution  Hamiltonian 


Based  on  an  analysis  of  the  necessary  conditions  and  engineering 
judgment  the  solution  obtained  appears  optimal  for  the  system  under 
consideration. 
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4.  Numerical  Considerations 

The  principal  problem  associated  with  CMG  systems  for  spacecraft 
attitude  control  is  the  well  known  geometric  singularity  conditions  This  refers  to  a 
condition  in  which  no  torque  is  generated  for  a  commanded  control  torque  in  a 
certain  “singular”  direction.  Mathematically,  this  occurs  when  the  A  matrix 
defined  in  equation  (6.5)  is  singular  or  rank  deficient.  By  solving  for  the  optimal 
control  vector  u*  to  minimize  an  objective  function  we  have  required  that  the 

control  solution  satisfy  equation  (6.5)  as  part  of  the  constraints  in  the  problem 
formulation.  As  a  result  the  possibility  of  the  matrix  A  being  singular  is 
eliminated.  The  condition  number  of  A  is  an  indication  of  how  close  this  matrix  is 
to  being  singular^.  Condition  number  is  defined  as  the  ratio  of  the  largest  to  the 
smallest  singular  value  of  the  matrix.  The  condition  number  of  an  identity  matrix 
is  one  and  the  condition  number  of  a  matrix  approaches  infinity  as  the  matrix 
becomes  singular.  The  condition  number  of  the  A  matrix  is  shown  in  Figure  88. 
It  is  clear  that  the  matrix  is  well  behaved  throughout  the  maneuver. 


Figure  88  Condition  Number  of  CMG  Control  Solution  Jacobian  Matrix 
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The  far-reaching  effects  of  this  singularity  free  solution  will  be  discussed  in  the 
next  chapter. 


5.  Conclusions 

In  this  section  the  open  loop  time-optimal  control  for  an  asymmetric 
spacecraft  equipped  with  four  CMGs  arranged  in  a  pyramid  configuration  was 
determined.  The  candidate  solution  was  determined  by  propagation  to  be  a 
feasible  solution  to  the  problem.  The  optimality  of  the  solution  was  validated 
through  an  analysis  of  the  Hamiltonian  minimization,  switching  functions  and  the 
behavior  of  the  Hamiltonian.  Therefore,  we  conclude  that  the  solution  is  feasible 
and  meets  the  necessary  conditions  for  optimality.  Additionally,  it  was  shown 
that  the  control  solution  is  completely  free  of  singularities.  This  leads  us  to 
conclude  that  if  the  control  solution  can  be  computed  in  real-time,  then  singularity 
avoidance  in  control  moment  gyro  actuator  based  systems  is  unnecessary. 
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VII.  CLOSED  LOOP  CONTROL 


A.  INTRODUCTION 

Classical  closed-loop  control  techniques  are,  in  general,  ill  suited  for  slew 
maneuver  time-optimal  control.  This  is  due  in  part  to  the  nonlinear  nature  of  the 
problem  and  the  inability  of  classical  control  techniques  to  generate  time-optimal 
solutions.  Additionally,  control  solutions  must  be  resolved  into  actuator  input 
parameters  to  generate  the  desired  torque.  This  method,  which  requires 
inverting  the  actuator  dynamics  equation,  has  led  to  singularity  problems  that 
continue  to  plague  engineers  as  evidenced  by  ongoing  researches, 4, 5. 

This  chapter  addresses  the  computation  of  real-time  optimal  controls  for 
spacecraft  slew  maneuvers.  Prior  research  efforts  were  focused  on  minimizing 
deviations  from  a  nominal  optimal  trajectory.  These  neighboring  optimal  control 
(NOC)  laws,  presented  by  Bryson  and  others6’7  have  been  used  in  a  variety  of 
applications. 

More  recently  a  two  degree-of-freedom,  non-linear  optimal  control  system 
architecture  was  suggested  by  Strizzi,  Yan,  et  al.8’9  This  concept  relied  on  the 
computational  power  of  pseudospectral  methods  to  develop  optimal  controls  and 
NOC  laws. 

In  this  work  optimal  control  solutions  to  the  non-linear  plant  are  generated 
to  implement  a  sampled-data  feedback  control  law.i°  Previous  optimal  solutions 
are  used  as  guess  for  subsequent  re-optimized  trajectories  and  significant 
reductions  in  computation  time  are  demonstrated. 

B.  CLASSICAL  CLOSED  LOOP  CONTROL 

A  basic  diagram  of  a  classical  closed  loop  spacecraft  attitude  control 
system  is  shown  in  Figure  89.  In  this  type  of  model  sensors  determine  the 
attitude  and  attitude  rates.  These  are  in  some  way  compared  with  desired 
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Motion 


Figure  89  Classic  Closed  Loop  Spacecraft  Attitude  Control  System 

attitude  and  rates  and  the  difference  or  error  signal  is  used  to  determine  a 
desired  control  torque  to  reduce  the  error  signal  to  zero.  Typical  approaches  for 
control  command  generation  involve  using  Euler  angle  errors  for  small  rotations 
and  either  a  direction  cosine  error  matrix  or  a  quaternion  error  vector  for  large 
angular  rotations.11  In  either  case  proportional  plus  derivative  gains  are  imposed 
in  order  to  achieve  desired  performance.  This  type  of  command  generation  has 
two  inherent  problems.  First,  the  command  generation  performs  no  better  than 
an  eigenaxis  maneuver  and  in  fact  will  seldom  perform  so  well.  Second, 
generating  a  desired  control  torque  for  use  as  an  actuator  input  has  led  to 
numerous  singularity  problems  of  which  the  best  known  is  the  control  moment 
gyro  problem. 


1 .  Error  Signal  Command  Generation 

In  previous  sections  we  have  seen  that  the  time-optimal  solution  is  a 
function  of  the  magnitude  of  the  maneuver  and  the  spacecraft  moment  of  inertia. 
While  error  signal  command  generation  is  capable  of  generating  the  eigenaxis 
path  the  command  signal  is  not  the  familiar  bang-bang  we  have  come  to  expect 
from  time-optimal  solutions.  As  the  spacecraft  approaches  the  desired  attitude 
the  command  generated  from  the  error  signal  approaches  zero.  Additionally,  the 
command  generation  algorithm  does  not  utilize  spacecraft  moment  of  inertia 
information  in  developing  the  command.  Modern  control  system  engineers  have 
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made  use  of  feedforward  signals  and  other  more  imaginative  techniques  to 
improve  slew  performance  however;  these  require  an  a  priori  knowledge  of  the 
maneuver  and  control  solution.  This  is  in  effect  a  step  back  from  the  autonomy 
that  we  seek  to  achieve. 

2.  Singularity  Issues 

We  have  seen  that  the  command  generation  algorithm  outputs  a  desired 
torque  output  from  the  actuator  based  on  the  error  signals  without  regard  to  the 
method  of  torque  generation.  This  desired  torque  is  sent  to  the  actuator  as  a 
commanded  torque  signal  to  be  accomplished.  The  question  to  consider  is 
under  what  circumstance  is  the  commanded  torque  actually  output  from  the 
actuator  and  what  are  the  effects  on  the  desired  maneuver  when  the 
commanded  torque  is  not  output  from  the  actuator. 

Consider  an  ideal,  unlimited  torque  thruster  which  is  mathematically 
modeled  as: 

"^"signal, 

”^signal2  —  ~^out2 
”^signal3  _  J~out$ 

In  this  case  the  input  signal  (Tsignal),  results  in  the  corresponding  output  torque 
{Tout).  Therefore,  given  the  desired  actuator  output  torque  {Tout ),  from  command 
generation,  we  can  solve  for  the  required  input  signal  to  the  actuator  (  Tsignal )  as, 

1  0  0 
0  1  0 
0  0  1 

The  actuator  dynamics,  represented  in  matrix  form,  in  this  case  as  the  identity 
matrix,  must  be  inverted  to  resolve  the  desired  output  torque.  In  this  case  the 
commanded  torque  generates  the  output  toque  desired.  Unfortunately,  actuator 
dynamics  are  seldom  this  simple.  As  a  result  we  define  an  actuator  singularity  as 
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a  condition  in  which  no  control  torque  is  generated  for  the  commanded  torque 
along  a  certain  direction.  Mathematically,  this  is  a  condition  in  which  the  actuator 
dynamics  matrix  is  singular. 

Consider  the  case  of  the  asymmetric  spacecraft  subject  to  CMG  control. 
This  case  was  examined  in  detail  in  Chapter  VI.  Since  we  have  determined  that 
the  error  command  generation  algorithm  is  unsuited  for  generating  bang -bang 
controls  suppose  that  we  elect  to  feed  forward  the  bang-bang  control  structure 
shown  in  Figure  90.  This  control  torque  on  the  spacecraft  would  result  in  a  135 
degree  slew  about  the  x-axis  for  NPSAT  1 . 


Figure  90  Asymmetric  Eigenaxis  Rotation  Torque  Control  Solution 

In  this  case  the  torque  command  signal  becomes  the  input  to  the  CMG  actuator 
in  our  spacecraft  model  (Figure  91 ).  The  input  or  commanded  torque  signal  must 
be  resolved  into  a  gimbal  angle  rate.  This  is  generally  accomplished  in 
accordance  with  a  command  pseudoinverse  steering  logic12.  In  this  steering 
logic  for  the  commanded  torque  input  (tv),  the  CMG  momentum  rate  command 
( h )  is  defined  as: 
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h  =  -U  -COX/7 

and  the  gimbal  rate  command  (8  )  is  then  obtained  by:* 

5  =A+h  =  AT(AATyh 


CMG  Open  Loop  Command  Propagation 


To  Workspace 


Figure  91  Spacecraft  Model  with  CMG  Actuators 


where  the  matrix  A  is  a  function  of  the  gimbal  angles  (5  )  and  was  previously 
defined  in  Chapter  VI  (Equation  (6.5)).  If  the  matrix  A  becomes  rank  deficient  for 
certain  sets  of  gimbal  angles  the  pseudoinverse  fails  to  exist  and  singular  states 
are  encountered.  Previously,  we  showed  that  the  time-optimal  CMG  solution 
resulted  in  well  behaved  condition  numbers  for  the  matrix  A  (See  Figure  88).  In 
this  case  the  bang-bang  control  solution,  when  used  as  the  desired  torque  on  the 
spacecraft,  produced  the  A  matrix  condition  numbers  shown  in  Figure  92. 
These  results  are  four-orders  of  magnitude  larger  than  those  obtained  in  the 


*  A+  is  also  called  the  Moore-Penrose  inverse  of  A  or  the  generalized  inverse  of  A. 
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Figure  92  Condition  Number  of  CMG  Jacobian  Matrix  with  Independent 

Torque  Solution 


Figure  93  Gimbal  Angle  Rates  From  Pseudoinverse  Steering  Logic 


previous  section.  As  a  result  of  the  poor  condition  of  this  matrix  the  algorithm 
attempts  to  drive  the  gimbal  angle  rates  to  infinity.  Gimbal  angle  rates  resulting 
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from  this  simulation  are  shown  in  Figure  93.  The  zoomed  in  view  of  the  gimbal 
angle  rates  from  4-6  seconds  is  shown  in  Figure  94. 


Figure  94  Gimbal  Angle  Rates  from  Pseudoinverse  Steering  Logic 

(Close-up) 

Gimbal  angle  rates  are  driven  towards  infinity  as  the  singular  states  are 
approached.  Clearly,  the  singularity  is  not  caused  by  the  bang -bang  structure  of 
the  input  torque.  Comparing  the  commanded  torque  signal  (Figure  90)  with  the 
condition  number  plot  (Figure  92)  or  the  gimbal  angle  rates  (Figure  93)  indicates 
that  the  singular  states  are  encountered  prior  to  the  control  switch. 

Singularity  avoidance  is  an  unintended  and  yet  significant  benefit  of  the 
time-optimal  control  algorithm.  In  effect  the  singularity  problem  is  avoided  both 
mathematically  and  physically  by  requiring  that  the  time-optimal  control  vector 
satisfy  the  actuator  constraints  in  the  problem  formulation  such  that, 

u*  =  -[>4(8  )]5  -oox/7(5) 

The  actuator  dynamics  are  solved  iteratively  from  left  to  right  vice  right  to  left. 
The  previous  matrix  inversion  is  eliminated  and  the  solution  is  optimal  with  regard 
to  the  actuator  capabilities. 
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In  this  example  we  have  examined  the  CMG  singularity  problem  since  it  is 
well  known.  However,  the  same  logic  could  be  applied  to  any  actuator  system 
that  is  not  an  ideal  torque  generating  device.  Magnetic  torque  rods  represent  a 
second  well  known  system  for  which  the  actuator  has  known  singularities.  In  this 
case,  a  desired  torque  output  from  the  torque  rod  must  be  resolved  in  the 
magnetic  dipole  moment  of  the  torque  rod.  The  torque  output  of  a  magnetic 
torque  rod  is  given  by: 


f  =  mxB 


'  0 

b3 

-b; 

>" 

f  = 

-B3 
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-B2 

-Si 
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Therefore,  given  the  desired  torque  output  we  have: 
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(7.2) 


The  inverse  of  equation  (7.1)  does  not  exist.  Therefore  it  is  impossible  to 
achieve  control  of  the  three  body  axes  by  means  of  the  created  magnetic  dipole 
vector  in  the  satellite13.  However,  using  the  time-optimal  control  formulation  as 
demonstrated  in  Chapter  III  and  Chapter  IV  we  require  that  the  optimal  control 
vector  (m)  satisfy  equation  (7.1).  Again,  the  singular  condition  is  avoided  by 
eliminating  the  matrix  inversion  and  iteratively  solving  the  problem  from  right-to- 
left. 


C.  BELLMAN  S  PRINCIPLE  OF  OPTIMALITY 

Bellman’s  principle  of  optimality  is  a  powerful  tool  when  applied  to  closed 
loop  time-optimal  control.  It  states  that  given  an  optimal  trajectory  from  a  point  A 
to  a  point  B,  the  trajectory  to  point  B  from  a  point  C  lying  on  the  optimal  trajectory 
is  also  optimal14. 
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Figure  95  Principle  of  Optimality 


In  this  work  this  principle  serves  three  purposes.  First,  it  verifies  optimality 
in  the  absence  of  costate  information.  Optimal  trajectories  are  obtained  and 
verified  by  recalculating  the  trajectory  from  an  intermediate  point  on  the 
trajectory.  Second,  in  the  presence  of  a  disturbance  torque,  this  principle  allows 
the  use  of  previous  optimal  solutions  as  guesses  for  subsequent  trajectories. 
Since  the  disturbed  state  is  near  the  optimal  state  using  the  previously  calculated 
state  trajectory  as  a  guess  for  the  re-optimized  state  trajectory  significantly 
reduces  the  computation  time.  Finally,  it  allows  us  to  evaluate  the  accuracy  of  a 
solution  based  on  a  finite  number  of  LGL  points. 

Shown  in  Figure  96  is  the  now  familiar  quaternion  trajectory  for  the 
asymmetric  spacecraft  time-optimal  slew  maneuver.  The  intermediate  point 
selected  is  on  the  optimal  trajectory  and  marked  with  a  bold  square.  The 
subsequent  solution  is  overlaid  with  “+”  marks. 


135 


Figure  96  Asymmetric  Time-optimal  Quaternion  Solution  -  Bellman’s 

Principle  of  Optimality 

The  subsequent  solution  clearly  obeys  Bellman’s  Principle  of  Optimality.  While 
this  may  not  be  sufficient  to  mathematically  demonstrate  the  optimality  of  the 
solution  it  can  certainly  be  considered  a  necessary  condition.  The  original 
solution  was  obtained  using  30  LGL  points  (nodes)  and  was  based  on  a  two- 
point  guess  defined  by  the  initial  and  final  conditions.  The  subsequent  solution 
also  used  30  nodes  but  was  provided  with  a  guess  which  consisted  of  15  points 
from  the  previous  solution.  The  computation  timest  as  determined  by  the 
MATLAB  tic  and  toe  commands  are  shown  in  Table  9 . 


t  All  computations  were  performed  on  a  Pentium  4  Windows  based  operating  system  at  3.06 
GHz  with  512  MB  RAM  using  MATLAB  6.5.0. 1  Release  13.  Emphasis  is  placed  on  the  relative 
vice  absolute  computational  time. 


136 


Computation  Time 

Original  Solution 

Subsequent  Solution 

Reduction,  % 

42.5310  seconds 

2.6250  seconds 

93.83 

Table  9  Bellman  Principle  Efl 

fect  on  Computational  Time 

The  accuracy  of  the  solution  obtained  is  related  to  the  number  of  nodes 
used  in  the  algorithm.  The  number  of  nodes  also  has  a  significant  effect  on 
computational  time.  Therefore,  we  strive  to  choose  a  number  of  nodes  that 
meets  our  accuracy  goal  without  creating  excessively  long  computational  times. 
As  the  number  of  nodes  increases,  the  values  returned,  states,  costates,  etc. 
converge  so  that  the  difference  between  solutions  approaches  zero.  Therefore,  if 
the  number  of  nodes  is  increased  and  the  behavior  of  the  solution  changes 
considerably  the  original  solution  is  shown  to  be  inaccurate.  In  Figure  97,  a  low 
node  (6  nodes)  solution  is  overlaid  with  a  high  node  solution  (30  nodes)  from  an 


137 


intermediate  point  on  the  optimal  trajectory.  The  significant  change  visible  in  the 
state  trajectory  and  final  time  are  indications  that  the  low  node  solution  has 
questionable  accuracy. 

D.  REAL-TIME  OPTIMAL  CONTROL 

In  this  section  we  provide  two  examples  of  closed  loop  time-optimal 
control.  In  each  case  the  simulation  is  as  shown  in  Figure  98.  The  primary 
tuning  parameter  for  spacecraft  performance  is  the  open-loop  propagation 
interval.  Ignoring  computational  delay,  the  state  is  sampled  and  an  open  loop 
control  solution  is  generated  time-indexed  to  the  predicted  spacecraft  state.  This 
solution  is  propagated  for  a  set  interval,  the  propagation  interval,  after  which  the 
state  is  re-sampled  and  a  re-optimized  trajectory  is  computed.  Currently, 
computation  times  are  recorded  for  comparison  and  evaluation  but  ignored  in  the 
simulation. 


Figure  98  Closed  Loop  Simulation  Model 


1 .  Asymmetric  Spacecraft  with  Idealized  Actuator  Control 

The  time-optimal  reorientation  of  an  asymmetric  spacecraft  with  idealized 
actuator  control  was  addressed  in  detail  in  Chapter  III.  The  spacecraft  and 
orbital  parameters  remain  unchanged  for  the  closed  loop  simulation  and  are 
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available  in  Chapter  III.  A  constant  disturbance  torque  has  been  added  to  the 
model.  This  disturbance  torque  will  perturb  the  model  off  the  optimal  trajectory 
and  require  re-optimization.  Since  we  know  from  the  previous  chapter  that  the 
time-optimal  open  loop  maneuver  is  approximately  6  seconds  the  disturbance 
torque  has  been  amplified  from  what  might  be  considered  normal  for  graphical 
visualization.  The  disturbance  torque  has  been  selected  as: 

Md  =0.108  N.m 

dx 

Md  =0.108  N*m 

ux 

Md  =0.031 7  N*m 

ux 

The  open  loop  time-optimal  control  solution  and  the  propagated  solution  in 
the  presence  of  this  disturbance  torque  are  shown  in  Figure  99.  The  solid  lines 
indicate  the  optimal  trajectory;  the  propagated  trajectory  is  marked  as  “+’.  It  is 


Quaternion  History  Solution  &  Propagation  Angular  Rate  History  Solution  &  Propagation 


Figure  99  Asymmetric  Spacecraft  Time-optimal  State  Solution  and 

Propagation 

clear  that  in  presence  of  a  disturbance  torque  the  open  loop  solution  would  not 
deliver  acceptable  performance. 

The  propagation  interval  for  this  closed  loop  simulation  was  selected  as  1 
second.  The  results  of  the  first  interval  of  propagation  are  shown  in  Figure  100 
and  Figure  101.  At  this  time  the  spacecraft  state  is  sampled  and  the  time-optimal 
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Figure  100  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  1 


Figure  101  Asymmetric  Spacecraft  Angular  Rate  Propagation  Interval  1 


solution  is  computed  from  the  current  attitude  and  rate  to  the  desired  attitude  and 

rate.  In  the  next  two  figures,  figure  102  and  Figure  103)  the  re-optimized 

trajectory  is  shown  along  with  the  original  trajectory.  The  basic  structure  of  the 

optimal  state  trajectories  is  the  same  and  the  end  point  conditions  are  met. 
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Figure  102  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  2 


Figure  103  Asymmetric  Spacecraft  Angular  Rate  Propagation  Interval  2 

The  sample  and  update  cycle  continues  from  Figure  104  through  Figure 
113.  Each  time  the  state  is  sampled  the  time-optimal  solution  is  generated  and 
propagation  continues  with  the  updated  solution. 
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Figure  104  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  3 


Figure  105  Asymmetric  Space  craft  Angular  Rate  Propagation  Interval  3 


With  each  update  it  is  clear  that  a  re-optimized  trajectory  is  generated  from  the 
current  state  of  the  spacecraft  to  the  final  state  which  remains  a  constant. 
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Figure  106  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  4 


Figure  107  Asymmetric  Spacecraft  Angular  Rate  Propagation  Interval  4 
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Quaternion  Solution  and  Propagation 


Figure  108  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  5 


Figure  109  Asymmetric  Spacecraft  Angular  Rate  Propagation  Interval  5 
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Figure  110  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  6 


Figure  111  Asymmetric  Spacecraft  Angular  Rate  Propagation  Interval  6 
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Figure  1 12  Asymmetric  Spacecraft  Quaternion  Propagation  Interval  7 


Figure  113  Asymmetric  Spacecraft  Angular  Rate  Propagation  Interval  7 


The  completed  maneuver  is  shown  in  Figure  112  and  Figure  113.  The 
time  required  to  complete  the  135  degree  reorientation  about  the  x-axis  in  the 
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presence  of  the  disturbance  torque  is  7.2652  seconds.  The  original  solution 
time,  modeled  without  a  disturbance  was  6.1825  seconds.  Each  subsequent 
solution  calculated  from  the  sampled  state  met  the  feasibility  and  optimality 
criteria  that  were  established  for  this  problem  formulation  in  Chapter  III. 

Computation  time  is  related  to  propagation  interval  in  that  longer 
propagation  intervals  in  the  disturbance  field  cause  the  spacecraft  to  deviate 
further  from  the  optimal  trajectory.  This  results  in  longer  computation  times.  The 
computation  times  for  the  simulation  are  shown  in  Figure  114.  The  average 
value  of  the  update  solutions  (solutions  2  through  7)  is  1 1 .9296  seconds.  This  is 
due  to  the  large  disturbance  torque  modeled  to  aid  in  the  visual  clarity  of  the 
graphics.  When  the  disturbance  torque  is  modeled  as  something  more  realistic 
for  the  spacecraft  and  orbital  parameters  the  previous  disturbances  are  reduced 
by  5  orders  of  magnitude.  The  resulting  update  computation  times  then 
decrease  to  an  average  value  of  4.0566  seconds. 


Figure  114  Asymmetric  Spacecraft  Idealized  Actuator  Time-optimal 

Solution  Computation  Time 
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a.  Numerical  Considerations  and  Notes 

Simulated  spacecraft  states  were  numerically  generated  by  a 
MATLAB  ODE45  propagation  subroutine.  The  generated  control  vector  plus  the 
disturbance  torque  were  used  as  input;  the  spacecraft  state  was  extracted 
following  propagation  to  simulate  a  sensor  reading.  In  some  cases  due  to 
numerical  approximations  and  truncation  errors  the  norm  of  the  quaternion  was 
not  equal  to  one  to  the  precision  required  for  the  optimization  algorithm  (DIDO)  to 
determine  a  feasible  solution.  In  order  to  resolve  this  problem  the  sampled-state 
quaternion  vector  was  normalized  as  follows: 


'nomal  |g||2 

Additionally,  Figure  114  indicates,  and  it  was  observed,  that  the 
computation  time  for  the  final  re-optimization  of  a  trajectory  was  in  general  longer 
than  previous  computation  times.  Since  this  algorithm  was  originally  conceived 
for  large  angle  slew  maneuvers  constant  scaling  parameters  were  applied  based 
on  the  maneuver  open  loop  solution.  As  the  maneuver  approaches  completion 
and  the  trajectory  is  re-optimized  the  magnitude  of  the  problem  has  changed  to 
the  extent  that  the  original  scaling  parameters  are  poorly  chosen.  The  problem 
requires  a  scaling  parameter  adjustment  based  on  the  remaining  horizon. 

2.  Asymmetric  Spacecraft  with  Magnetic  Torque  Rod  Control 

The  time-optimal  reorientation  of  an  asymmetric  spacecraft  under 
magnetic  torque  control  was  treated  in  detail  in  Chapter  IV.  Once  again  the 
spacecraft  and  orbital  parameters  are  based  on  the  NPSAT  1  small  satellite 
program  and  found  in  Table  5. 

The  closed  loop  simulation,  while  made  more  complex  by  the  time-varying 
magnetic  field,  is  accomplished  as  shown  in  Figure  98.  The  constant  disturbance 
torque  selected  for  propagation  is: 
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Md  =5.4x1  O'6  N.m 

ux 

Md  =5.508x10  6  N*m 

y 

Md  =6.34x10  7  N.m 

dz 

The  larger  disturbance  torque,  used  in  the  previous  example  for  visual  clarity, 
would  overpower  the  torque  available  from  the  system  and  make  the  spacecraft 
uncontrollable.  Recall  that  the  time-optimal  open-loop  maneuver  required 
approximately  273  seconds.  The  state  trajectories  resulting  from  the  propagation 
of  the  open  loop  solution  in  the  presence  of  the  disturbance  torque  are  shown  in 
Figure  115  and  Figure  116.  The  end  point  error,  although  difficult  to  see  in  the 
graphical  presentation  is  significant  enough  to  warrant  correction.  A  comparison 
of  the  actual  vs.  desired  state  is  shown  in  Table  1 0. 


Figure  1 15  Magnetic  Torque  Open  Loop  Solution  Propagated  with 

Disturbance  Torque 


149 


Figure  116  Magnetic  Torque  Open  Loop  Solution  Propagated  with 

Disturbance  Torque  (2) 


Propagation  Error  with  Disturbance 

State 

Desired 

Actual 

Error 

ql 

0.9239 

0.9216 

0.0023 

q2 

0 

-0.01 

0.01 

q3 

0 

0.0206 

0.0206 

q4 

0.3827 

0.3879 

0.0052 

Omega  1 

0.00 

-1.00E-04 

1.00E-04 

Omega  2 

7.72E-04 

7.00E-04 

7.25E-05 

Omega  3 

7.72E-04 

1.30E-03 

5.28E-04 

Table  1 0  State  Trajectory  Propagation  Error  with  Disturbance 


A  propagation  interval  of  15  seconds  was  selected  for  this  simulation 
based  on  the  open  loop  maneuver  time  and  the  magnitude  of  the  disturbance 
torque.  As  before,  following  propagation  the  spacecraft  state  is  sampled  and  the 
trajectory  is  re-optimized  from  the  current  state. 
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DIDO  has  internal  constraints  which  define  the  accuracy  with  which  the 
optimal  trajectory  must  reach  the  defined  end  points.  In  order  to  facilitate  the 
solution,  box  constraints  were  applied  to  the  desired  end  state.  These  form  an 
upper  and  lower  bound  on  the  end  state  such  that  and  acceptable  end  state  is  ± 
epsilon  from  the  exact  end  state.  For  this  simulation  epsilon  was  set  as, 

q£  =0.006 
coE  =0.0001 

The  end  state  box  constraints  appeared  to  increase  computational  speed. 
In  Figure  117,  the  first  update  to  the  disturbed  quaternion  trajectory  is  shown  and 
appears  to  complete  the  maneuver  faster  than  the  original  undisturbed  solution. 
This  is  not  surprising.  The  open  loop  solution  was  formulated  to  the  end  point 
and  the  subsequent  solution  was  formulated  to  the  end  point  with  the  addition  of 
the  box  constraint.  Additionally,  it  is  possible  that  the  constant  disturbance 
torque  is  initially  providing  some  benefit  to  maneuver  speed.  Flowever,  the  final 
closed  loop  state  trajectories  figure  118  and  Figure  119)  demonstrate  that  to 
meet  the  established  accuracy  requirements  in  the  presence  of  the  disturbance 
torque  requires  additional  time. 
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Quaternion  Solution  and  Propagation 


Figure  117  Magnetic  Torque  Quaternion  Propagation  Interval  2 


Figure  118  Magnetic  Torque  Quaternion  Propagation  -  Maneuver 

Completion 
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Figure  119  Magnetic  Torque  Angular  Rate  Propagation  -  Maneuver 

Completion 

A  comparison  of  the  open  loop  maneuver  with  the  disturbed  closed  loop 
maneuver  shows  that  the  disturbance  increases  the  overall  maneuver  time. 

The  disturbance  free  simulation  required  273.4117  seconds  to  complete  the  135 
degree  x-axis  slew.  The  close  loop  maneuver  in  the  presence  of  the  disturbance 
torque  required  285.3174  seconds.  The  overall  accuracy  improvement  is  shown 
in  Table  11.  The  improvement  appears  significant  enough  to  warrant  the 
application  of  closed  loop  control.  Again,  all  subsequent  solution  updates  met 
the  feasibility  and  optimality  necessary  conditions  that  were  established  in 
Chapter  IV. 


Closed  Loop  Propagation  Error 

State 

Desired 

Actual 

Error 

Error  Reduction,  % 

ql 

0.9239 

0.9263 

0.0024 

4.35 

q2 

0 

0.0042 

0.0042 

58.00 

q3 

0 

0.006 

0.006 

70.87 

q4 

0.3827 

0.3767 

0.006 

15.38 

Omega  1 

0.00 

-1.00E-04 

1.00E-04 

0.00 

Omega  2 

7.72E-04 

8.00E-04 

2.75E-05 

62.01 

Omega  3 

7.72E-04 

7.00E-04 

7.25E-05 

86.26 

Table  1 1  Closed  Loop  State  Propagation  Results 
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As  before,  computation  time  is  the  key  to  establishing  the  feasibility  of 
closed  loop  real-time  optimal  control.  The  computation  of  magnetic  torque  time- 
optimal  maneuvers  is  significantly  complicated  by  the  time-varying  magnetic 
fields.  This  is  illustrated  by  the  long  computation  time  required  for  the  open  loop 
solution  shown  in  Figure  120.  The  open  loop  solution  was  based  on  a  two-point 
guess,  the  end  points,  and  calculated  over  30  nodes.  The  subsequent  solutions 


Figure  120  Magnetic  Torque  Time-optimal  Maneuver  Computation  Time 

were  also  calculated  over  30  nodes  but  used  the  previous  state,  control  solution 
as  a  guess.  The  time  required  for  the  open  loop  solution  was  324.9680  seconds, 
subsequent  solutions  had  an  average  time  of  12.3407  seconds.  This  represents 
a  computation  time  reduction  of  96.2  percent. 


E.  CONCLUSIONS 

This  section  has  demonstrated  that  classical  feedback  control  theory  is  ill 
equipped  to  handle  time-optimal  maneuvers.  More  recent  techniques  employing 
neighboring  optimal  control  laws  have  improved  slew  performance  but  still 
require  inverting  the  actuator  dynamics  and  therefore  require  singularity 
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avoidance.  By  requiring  that  the  control  solution  satisfy  the  actuator  dynamic 
constraint  equations  singularity  avoidance  is  automatically  eliminated. 

Strizzi  et  al. 15  noted,  and  engineering  common  sense  agrees,  that  if  open 
loop  control  solutions  can  be  generated  in  real-time  they  become  equivalent  to 
feedback  control  laws.  Computational  power  and  numerical  methods  have  come 
together  to  make  this  a  reality. 

The  two  numerical  examples  presented  in  this  section  illustrate  the 
concept  of  fast  and  slow  dynamics.  In  the  idealized  actuator  example,  a 
computational  time  of  10  seconds  would  not  suffice  for  real-time  control. 
However,  that  computation  time  applied  to  the  magnetic  torque  example  would 
seem  an  excessive  requirement.  The  computational  speed  required  to  achieve 
real-time  optimal  control  is  a  product  of  both  the  plant  dynamics  and  desired 
performance  characteristics. 

Empirically  it  was  observed  that  the  time  required  to  calculate  the 
subsequent  trajectories  is  a  function  of  two  factors,  deviation  from  optimal 
trajectory  and  amount  maneuver  remaining.  As  the  frequency  at  which  updates 
are  made  increases,  the  amount  of  deviation  from  the  optimal  trajectory  due  to 
disturbance  torques  and  plant  uncertainties  necessarily  decreases.  Therefore, 
computation  time  required  decreases  as  update  frequency  increases.  Second, 
as  the  maneuver  approaches  completion  the  computation  time  required 
asymptotically  approaches  a  minimum  value. 

Using  previous  solutions  to  jump  start  the  updated  solution  was  shown  to 
significantly  reduce  solution  computation  time.  Further  reductions  in  actual 
computation  time  are  predicted  by  reducing  operating  system  overhead  and 
employing  field  programmable  gate  arrays. 
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VIII.  FUTURE  WORK 


The  work  in  this  thesis  is  by  no  means  ready  for  on-orbit  testing.  What 
follows  are  areas  of  improvement  and  future  research  that  time  precluded. 


A.  ALGORITHM  IMPROVEMENT  AND  OPTIMIZATION 

This  algorithm  was  originally  conceived  for  large  angle  slew  maneuvers. 
To  that  end,  numerical  scaling  was  used  to  adjust  all  parameters  to  within  the 
same  order  of  magnitude.  This  improved  the  accuracy  of  results  and  reduced 
computational  time.  However,  as  the  spacecraft  moment  of  inertia  was  varied 
there  were  conditions  under  which  separate  scaling  for  the  angular  rates  of  the 
body  axes  would  have  improved  performance.  A  scaling  format  based  on 
structure  array  variables1  would  simplify  coding  and  allow  separate  scaling  for 
each  variable. 

Additionally,  the  time  scaling  was  designed  for  the  initial  open  loop 
solution  and  held  constant  throughout  the  closed  loop  simulation.  Therefore,  as 
the  time-to-go  horizon  decreased,  the  problem  was  poorly  scaled  in  some  cases. 
As  the  spacecraft  state  is  sampled  for  each  re-optimization  the  time  scale  should 
be  adjusted  so  as  to  ensure  that  subsequent  solutions  are  well  scaled  for 
accuracy  and  to  reduce  computational  time. 

Though  great  emphasis  has  been  placed  on  computational  time  in  the 
realization  of  real-time  optimal  control,  the  algorithms  used  in  this  work  were  not 
designed  to  be  time-efficient.  Some  reconfiguration  of  the  loops  and  control 
structure  of  the  algorithms  would  undoubtedly  result  in  additional  computational 
speed. 

Finally,  though  not  the  subject  of  this  work,  increasing  the  speed  of  the 
reusable  software  package  DIDO  is  being  investigated.  Through  encoding 
analytic  Jacobian  significant  computational  speed  enhancements  are  possible.2-3 
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B.  COMPUTATIONAL  DELAY  MODELING 


Key  to  the  implementation  of  this  algorithm  is  the  computational  delay 
associated  with  the  solution  to  the  optimal  control  problem.  As  shown  in  Figure 
121  the  computational  delay  (Sc)  is  the  time  from  when  the  state  is  observed  at  ti 

until  the  time  when  the  re-optimized  control  solution  is  available  at  t2.  During  that 
time  the  state  is  continuing  to  propagate  under  the  previous  control  solution. 


Figure  121  Conceptual  Computational  Delay  Model 
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Bellman’s  principle  of  optimality  says  that  in  the  absence  of  a  disturbance 
torque  the  sampled  state  at  ti  would  lie  on  the  optimal  trajectory  and  the  updated 
control  solution  would  be  identical  to  the  original  solution.  However,  in  the 
presence  of  a  disturbance,  the  sampled  state  could  be  off  the  optimal  trajectory 
and  the  updated  control  solution  would  not  be  applied  until  some  non-zero  time 
later.  In  this  case,  the  re-optimized  solution  would  be  applied  to  a  system  whose 
state  has  propagated  to  another  point,  presumably  off  the  re-optimized  trajectory 
thus  inducing  some  error  into  the  system.  In  this  work,  the  computational  delay 
was  taken  to  be  negligibly  small.  In  fact,  the  effect  of  the  computational  delay  is 
a  function  of  both  the  spacecraft  dynamics  and  the  desired  performance. 
Numerical  modeling  and  simulation  of  the  effect  of  computational  delay  is 
necessary  to  determine  both  a  sample  rate  and  maximum  acceptable 
computational  delay. 

C.  EXTENSION  TO  HARDWARE 

Finally,  much  research  is  needed,  and  has  begun,  to  implement  these 
algorithms  in  hardware.  Implementation  via  field  programmable  gate  arrays 
(FPGA)  is  in  ongoing.  Through  the  efforts  of  others  this  project  will  see  air¬ 
bearing  testing  and  eventual  on-orbit  testing. 
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